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INTRODUCTION 

Progress  on  the  development  of  the  wave  drag  codes  and 
improvements  in  the  radiation  code  were  considerably  slowed  by 
the  changeover  to  the  UCLA  360/91  computing  system.  Our  un¬ 
familiarity  with  the  new  system,  inadequate  consulting  services, 
and  frequent  shutdowns  of  our  host  (the  UCSD  computing  center) 
contributed  to  long  delays  in  our  program  plan.  Most  of  the 
work  which  remains  to  be  completed  in  this  contract  would  have 
been  accomplished  during  this  report  period  f  we  had  been  oper¬ 
ating  on  our  local  computing  system. 

At  the  close  of  the  last  contract  period,  the  linear 
steady-state  model  had  been  formulated.  The  stress  integral 
had  been  separated  into  an  atmospheric  response  function  and 
a  topography  spectrum  function.  All  of  the  information  about 
the  atmosphere  had  been  contained  in  the  atmospheric  response 
function, and  all  the  information  about  the  topography  was  con¬ 
tained  in  the  topography  spectrum  function. 

Our  objective  in  this  contract  period  has  been  to  param¬ 
eterize  the  drag  induced  by  the  mountain  lee  waves  by  characteriz¬ 
ing  the  atmospheric  response  function  in  terms  of  the  atmospheric 
variables  provided  by  the  Rand  Global  Circulation  Model  (GCM) . 
Specifically,  these  variables  include  the  surface  pressure  in 
addition  to  the  temperatures  and  meridional  and  zonal  winds  at 
the  levels  of  the  GCM.  The  parameterization  of  the  wave  drag 
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in  forms  of  the  variables  from  the  two-level  model  required  as¬ 
sumptions  to  be  made  about  the  atmospheric  wind  and  temperature 
profiles.  Once  these  assumptions  were  made,  the  atmospheric  re¬ 
sponse  function  and  the  wave  drag  could  be  characterized  by  study¬ 
ing  the  response  of  many  atmospheres  taken  from  a  previous  run 
of  th1*  Global  Circulation  Model. 

The  program  during  the  past  six  months  involved  a  series  of 
tasks  all  leading  to  the  parameterization  of  the  lee  wave  drag.  A 
subroutine  will  be  delivered  to  the  Rand  Corporation  by  the  end  of 
the  contract  period  to  be  tested  and  subsequently  incorporated  into 

the  GCM.  The  below  tasks  were  completed  during  this  contract 
period: 

1.  We  became  familiar  with  the  ARPA  network  an:  made 
our  codes  compatible  with  the  system. 

2.  A  prescription  for  the  atmospheric  wind  and  tempera¬ 
ture  profiles  was  selected. 

3.  The  atmospheric  response  function  code  was  rewritten 
and  optimized. 

4.  The  code  for  calculating  the  stress  was  written. 

5.  The  topography  spectrum  function  was  run  and  studied 
for  varied  topographies. 

6.  The  atmospheric  response  function  was  characterized 
in  terms  of  data  from  the  two-level  model. 

7.  The  stress  integral  was  studied  as  a  function  of 
topography. 

8.  The  component  of  the  stress  in  the  direction  of  the 
meridional  winds  was  characterized,  and  a  subroutine 
based  on  an  analytic  model  has  been  formulated. 

The  same  problems  which  plagued  the^esoscale  program  also 
effected  the  work  on  the  radiation  codes.  Most  of  the  work  per¬ 
formed  during  this  period  involved  making  the  codes  compatible 
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with  the  IBM  3(0/91  system  and  performing  trial  calculations  to 
chock  out  these  codes.  Modifications  to  the  codes  were  made  in 
preparation  for  the  delivery  of  these  codes  to  the  Rand  Climate 
Dynamics  Group. 

1.  CHOICE  OF  WIND  AND  TEMPERATURE  PROFILES  IN  THE  ATMOSPHERE 

The  amplitude  of  the  waves  generated  on  the  leeward  side 
of  mountains  and  the  mechanism  for  the  dissipation  of  these 
waves  will  depend  strongly  on  the  atmospheric  wind  and  tempera¬ 
ture  profiles.  In  the  linearized  theory,  the  propagation  of 
wave  energy  in  the  wave  with  horizontal  wavenumbers  (k,l)  is  de¬ 
pendent  on  the  Scorer  parameter  profile  in  the  atmosphere.  The 
Scorer  parameter  l2(z)  =  N2/Un2  -  d2Un/dz2  is  a  function  of  the 
atmospheric  stability,  N,  and  the  velocity  component,  Un,  parallel 
to  the  wavevector.  The  amplitude  of  the  waves  that  the  atmosphere 
can  support  is  determined  by  a  solution  of  the  equation 

(d2w/dz2)  +  (l2(z)  -  k2)  w  =  0  (1.1) 

subject  to  the  boundary  conditions  imposed  by  the  model  used. 

2  2  2 

In  equation  (1.1),  k  =  k  +  1  .  The  local  magnitude  of  the  Scorer 
parameter,  1  (z) ,  with  respect  to  the  square  of  the  magnitude  of 
the  wavevector,  k  ,  will  determine  the  amplification  of  a  parti¬ 
cular  wavenumber.  When  the  wind  speed  in  the  direction  of  the 
wavevector  is  very  high,  horizontal  energy  and  momentum  is  pumped 
into  that  wavenumber,  and  the  wave  will  grow  in  amplitude  in  that 
layer.  Conditions  producing  amplification  of  wave  energy  (<2>l2(z)) 
also  occur  when  the  lapse  rate  approaches  the  adiabatic  lapse  rate. 
The  longer  waves  will  not  be  as  greatly  influenced  by  r\anges  in 
the  resonance  characteristics  of  the  lower  atmosphere,  and  much 
of  this  wave  energy  will  pass  through  the  troposphere  forcing  the 
natural  modes  of  the  stratosphere. 
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The  generation  and  dissipation  of  lee  waves  in  an  atmosphere 
wi]]  depend  upon  the  fine  structure  of  the  wind  field  and  the 
temperature  profile.  The  response  characteristics  of  the  atmos¬ 
phere  to  a  disturbance  with  wavenumber  (k,l)  is  determined  by 
the  local  stability  of  the  atmosphere,  the  wavelength  of  the 
disturbance,  and  the  angle  the  wavevector  makes  with  the  meri¬ 
dional  and  zonal  winds.  The  most  favorable  conditions  for  the 
generation  of  high  amplitude  waves  occur  with  strong  low-level 
winds  and  a  low-level  inversion.  The  dissipation  of  the  waves 
will  occur  at  regions  of  low  Richardson's  number,  that  is,  re¬ 
gions  of  high  stress.  The  effect  of  horizontal  stress  on  the 
wave  is  to  turn  the  wavevector  in  the  horizontal  direction.  (6) 

As  the  intrinsic  frequency  of  the  wave,  ui  —  U  < ,  goes  to  zero, 
the  amplitude  decreases  and  the  momentum  is  given  to  the  flow- 
field.  Hence,  an  accurate  description  of  the  momentum  trans- 
Port  to  the  atmosphere  requires  a  detailed  atmospheric  tempera¬ 
ture  ard  wind  velocity  profile. 

The  modeling  of  the  wavedrag  associated  with  mountains  is 
restricted  by  the  lack  of  vertical  and  horizontal  resolution  of 
a  global  circulation  model.  We  can  resolve  the  topography  to  the 
5  or  30  grid  data  available,  but  we  rely  on  the  atmospheric 
/ ariables  from  the  GCM  for  a  definition  of  the  atmospheric  pro¬ 
files  of  temperature  and  wind  velocity.  The  wind  vectors  and 
temperatures  are  resolved  to  only  the  4°  by  5°  grids  used  in  the 
GCM.  The  strong  low-level  inversions  (which,  in  any  case,  could 
not  be  resolved  in  the  vertical  by  a  two— level  model)  are  gener¬ 
ally  the  result  of  the  movement  of  a  warm  air  mass  over  a  cold 
air  mass.  The  resolution  of  frontal  activity  cannot  be  obtained 
explicitly  from  the  GCM.  However,  arguments  for  the  occurrence 
of  inversions  and  a  parameterization  of  the  inversion  level  may 
be  possible  from  a  study  of  the  movement  of  synoptic  scale  air 
masses  in  the  GCM.  We  ha^e  not  had  the  opportunity  to  attempt 
such  a  parameterization  in  this  contract  period,  but  an  approach 
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of  this  type  could  improve  the  horizontal  and  vertical  resolu¬ 
tion  of  the  atmospheric  conditions.  Our  lee-wave  model  is 
limited  by  the  resolution  of  the  topography  data  available  and 
not  by  the  4°  by  5°  grids  of  the  GCM.  Any  arguments  which 
would  improve  the  resolution  of  the  wind  and  temperature  data 
would  improve  the  treatment  of  the  wave  drag. 

The  wave  drag  model  that  we  are  applying  is  the  linear, 
steady-state,  adiabatic,  Boussinesq  formulation  of  Bretherton.  ^ 
The  momentum  flux,  in  this  model,  is  independent  of  height  in 
the  atmosphere  except  at  critical  layers  where  dissipation  of 
the  momentum  occurs.  The  boundary  conditions  of  the  model  do 
not  permit  reflection  of  waves  and  complete  absorption  of  a 
wave  will  o:cur  at  a  critical  layer.  The  response  of  the  at¬ 
mosphere  is  determined  by  a  solution  of  equation  (1.1)  subject 
to  the  boundary  conditions, 

3“  =  -/^2  -  i2(hiw 

dz 

or  ^  =  i  \/T2(h)  "  sgn  (Un)  w , 

applied  at  the  level  z  =  H  where  the  calculation  is  begun.  The 
first  condition  in  equation  (1.2)  is  applied  when  the  wavenumber 
k2  is  greater  than  the  Scorer  parameter  at  z  =  H  and  represents 
exponential  decay  of  the  wave  amplitude  above  the  cutoff  height. 
The  second  condition  is  applied  when  the  Scorer  parameter  is 
greater  than  the  wavenumber  and  represents  a  radiation  condi¬ 
tion  in  which  the  waves  continue  to  propagate  upward  through 

the  layer  at  z  =  H.  A  detailed  discussion  of  this  model  has 

(7  8 ) 

been  given  in  previous  reports.  ' 

The  effort  during  this  phase  of  the  contract  period  has 
been  directed  at  characterizing  the  stress  integral  and  deter¬ 
mining  to  what  degree  the  integrands  of  the  two  components  of 
the  stress  could  be  simplified.  In  characterizing  the  stress 
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in tuyral ,  we  have  chosen  to  fit  the  GCM  wind  speed  and  teirpcra- 
ture  data  at  the  two  level*-  with  prescribed  vertical  profiles 
which  will  be  described  subsequently.  The  Scorer  parameter  pro¬ 
file  for  each  wavevector  then  followed  from  the  definition  of 
the  Scorer  parameter  and  the  atmospheric  response  was  determined 
by  integration  of  equation  (1.1)  through  the  atmosphere  subject 
to  the  boundary  conditions  at  the  top  of  the  atmosphere  given 
in  equation  (1.2).  The  wind  profiles,  the  temperature  profiles, 
and  the  Scorer  parameter,  are  of  necessity  considerably  smoothed 
due  to  the  poor  resolution  of  the  GCM.  Since  the  Scorer  param¬ 
eter  is  considerably  smoothed,  we  might  have  alternately  chosen 
to  express  the  Scorer  parameter  in  term’s  of  the  mathematical 
forms  which  give  analytic  solutions  to  the  Scorer  equation.  Ap¬ 
plication  of  these  analytic  solutions  to  multilayer  models  for 
the  one-dimensional  problem  has  >  en  given  by  Palm  and  Foldvik 
and  Vergeiner. ^,10)  The  same  solutions  are  applicable  to  our 
problem  if  the  Scorer  parameter  is  treated  as  a  function  of  the 
wavevector  of  each  of  the  components  of  the  topography  spectrum 
function.  The  Scorer  equation  admits  analytic  solutions  when 
the  Scorer  parameter  is  treated  as  a  constant  through  the  layer, 
as  an  exponential  function  of  height,  or  if  it  .< s  taken  to  be 


l2(z) 


(z+a) 


(1.3) 


where  1Q  and  a  are  constants  determined  for  each  level.  The 
form  given  in  equation  (1.3)  has  been  applied  by  Vergeiner.  We 
will  discuss  in  Section  3  how  we  intend  to  use  the  analytic  solu¬ 
tions  in  the  Rand  GCM  subroutine  to  calculate  the  two  components 
of  the  wave  drag.  We  did  not  feel  justified  in  using  these  ana¬ 
lytic  models  during  this  phase  of  the  contract  period  because  of 
the  need  to  study  the  sensitivity  of  the  atmospheric  response 
function  to  the  changes  in  the  atmospheric  variables. 
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1 . 1  Mhe  Choice  of  the  Wind  Profile 

We  have  described  the  difficulties  involved  in  obtaining 
the  meridional  and  zonal  wind  profiles.  We  inspected  the  wind 
data  chat  were  available  on  National  Weather  Service  tapes, 
obtained  aata  available  on  microfilm  from  the  National  Center 
for  Atmospheric  Research,  and  reviewed  the  seasonally  averaged 

( 12 ) 

wind  profiles  available  from  the  National  Weather  Records  Center. 

We  also  consulted  with  meteorologists  to  obtain  suggestions  on 
how  the  wind  profiles  might  be  deduced  from  the  two-level  model. 

In  Figures  1  and  2  we  have  reproduced  typical  seasonally 
averaged  data  from  the  National  Weather  Records  Center  data  for 
the  midlatitudes  over  mountainous  terrain.  Each  component  is 
characterized  by  peak  speeds  near  the  tropopause  and  can  be  ap¬ 
proximated  by  a  second-degree  polynomial. 

Wind  speed  and  direction  data  obtained  from  the  National 
Center  for  Atmospheric  Research  are  reproduced  in  Figure  3. 

These  data  represent  soundings  over  the  Rocky  Mountains  taken 
over  a  period  of  approximately  one  hour. 

Finally,  in  Figures  4  to  6  we  have  reproduced  data  taken 
from  the  National  Weather  Service  tapes  for  the  meridional  and 
zonal  winds  and  the  temperature  profiles  over  West  Virginia. 

The  data  are  represented  by  the  dots  and  the  data  points  have 
been  fitted  with  a  smooth  curve. 

In  this  parameterization  study  the  solution  of  the  Scorer 
equation  was  obtained  by  an  integration  of  the  equation  through 
the  atmosphere  in  the  z-coordinate  system.  To  use  the  GCM  data 
required  a  transformation  from  the  c- coordinate  system  to  the 
z-coordinate  system  and  extrapolation  of  the  data  to  other  points 
in  the  atmosphere.  The  temperature,  given  at  the  o  =  1/4  and 
o  =  3/4  levels,  was  extrapolated  in  a  pressure  coordinate  system 
by  letting 


7 


bSS-R-74-2331 


-10  -3  -6  -4  -2  0  2  4  6  8  ] 

Meridional  fie  an  Uind  from  South  (Knots) 

•  Dec  -  Jan  -  Feb 
a  Mar  -  Apr  -  Aay 
d  Jun  -  Jul  -  Aug 
■  Sep  -  Oct  -  Nov 

Figure  2.  Seasonally  averaged  meridional  mean  wind  for  110°W 
45°N  from  U.S.  Dept,  of  Commerce  Technical  Paper  41 
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/  GCM  prescription 


Temperature  (°k) 


The  temperature  profile  for  the  West  Virginia  data 
taken  from  the  National  Weather  Service  tapes. 
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belov  the  200-mb  level  and  assuming  an  isothermal  layer  above 
this  level.  In  equation  (1.4),  =  0.286,  and  subscripts  1 

and  3  denote  the  variables  at  the  o  =  1/4  and  o  =  3/4  levels, 
respectively.  The  transformation  from  the  p-coordinates  to  the 
z-coordinates  followed  by  successive  application  of  the  hydro¬ 
static  equation  written  in  finite  difference  form.  The  pressure 
profile  obtained  in  this  manner  from  the  West  Virginia  data  is 
given  in  Figure  7. 

The  prefiles  of  the  meridional  aid  zonal  winds  were  ob¬ 
tained  by  a  parabolic  fit,  in  z-coordinates,  of  the  wind  data 
points  at  o  =  1/4,  o  =  3/4,  and  a  third  point  obtained  by  a  re¬ 
flection  of  the  o  =  1/4  data  point  about  the  200-mb  pressure 
level.  This  profile  closely  approximates  the  profiles  of  the 
seasonally  averaged  meridional  and  zonal  winds  in  the  midlati¬ 
tudes.  In  Figures  4  to  6  we  have  applied  the  extrapolation 
schemes  to  the  West  Virginia  wind  and  temperat  ire  profiles. 

In  this  particular  example,  very  high  surface  winds  are  pre¬ 
dicted  by  this  model.  We  did  not  correct  for  the  high  surface 
winds  because  in  general  this  viar  not  a  problem.  The  only  ad¬ 
vantage  of  this  prescription  for  the  winds  is  that  it  provides 
an  extrapolation  above  200  mb  which  is  not  given  by  the  linear 
extrapolation  in  the  GCM.  A  prescription  for  the  wind  field 
which  more  accurately  treats  the  variability  of  tl.e  winds  with 
height  would  require  a  lengthy  statistical  treatment  of  the 
available  wind  data  in  terms  of  the  two  levels  of  the  GCM.  The 
additional  information  gained  by  such  a  prescription  would  un¬ 
doubtedly  be  lost  in  the  parameterization  required  to  reduce 
the  computing  time  to  levels  compatible  with  the  GCM.  A  simple 
modeling  of  the  wind  field  is  consistent  with  the  requirements 
of  a  simple  mathematical  expression  for  the  Scorer  parameter  to 
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Figure  7.  The  pressure  profile  for  West  Virginia  from  the 
National  Weather  Service  tapes. 
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obtain  an  analytic  solution  of  the  Scorer  equation.  In  general, 
in  the  absence  of  very  stable  conditions,  the  Scorer  parameter 
will  be  much  larger  than  the  wavenumbers  resolvable  in  the  topo¬ 
graphy  spectrum.  Any  local  decrease  of  the  Scorer  parameter 
will  have  a  greater  influence  on  the  shorter  wavelengths  than 
on  the  long  wavelengths.  Other  investigators  have  found  that 
these  trapped  waves  do  not  increase  the  wave  drag  significantly . ^ 

1. 2  Choice  of  the  Temperature  Profile  and  Its  Effect  on  the 
Atmospheric  Response 

The  temperatures  at  the  a  =  1/4  and  a  =  3/4  levels  have 
been  extrapolated  to  other  levels  by  assuming  that  the  potential 
temperature  is  linear  in  pkl,  where  =  1  -  1/y  and  y  is  the 
ratio  of  specific  heats  for  the  atmosphere.  We  have  extended 
this  to  the  200-mb  level,  and  then  we  have  assumed  an  isothermal 
atmosphere  above  this  level.  The  extrapolation  of  the  tempera¬ 
ture  is  the  rame  as  used  in  the  GCM. 

The  temperature  profile  enters  explicitly  in  the  Scorer 
equation  through  the  stability.  It  also  enters  into  the  wind 
profile  through  the  transformation  of  the  a-cocrdinate  system 
to  the  z-coordinate  system  by  the  hydrostatic  equation.  The 
smoothing  of  the  temperature  profile  should  also  have  a  more 
significant  effect  on  the  shorter  wavelengths  than  on  the  long 
wavelengths . 

2.  BEHAVIOR  OF  THE  STRESS  INTEGRALS 

We  have  described  the  linear,  steady-state  Bretherton 

(78) 

model  in  detail  in  previous  reports'  '  and  have  derived  the 
stress  integrals  in  terms  of  the  atmospheric  response  function 
and  the  topography  spectrum  function.  In  this  section  we  de¬ 
scribe  the  results  of  parameterization  studies  intended  to  lead 
to  a  subroutine  for  the  mountain  lee  wave  drag  compatible  with 
the  Rand  Global  Circulation  Model.  We  have  considered  the  ten 


► 


SSS-R-74-2331 


topographic  grids  and  the  15  atmospheres  listed  in  Tables  1  and 
2.  All  of  the  grids  listed  in  Table  1  were  2°40‘  on  a  side. 

In  addition,  another  computer  run  was  made  of  the  stress  inte¬ 
gral  on  a  4°  by  5°  grid  with  the  northwest  corner  located  at 
12].°  w  37°  N.  The  smaller  grids  were  selected  for  the  economy 
of  calculating  fewer  values  of  the  atmospheric  response  func¬ 
tion.  The  atmospheres  were  taken  from  selected  mountainous 
areas  of  a  GCM  run.  Four  atmospheres  were  taken  from  this  GCM 
run  and  the  remaining  11  atmospheres  are  perturbations  or  per¬ 
mutations  of  these  four  sets  of  data. 

The  topography  spectrum  function,  the  atmospheric  re¬ 
sponse  function,  and  the  stress  integral  are  discussed  in  the 
next  three  sections.  Our  objective  was  to  parameterize  the 
stress  integral;  therefore,  much  of  the  discussion  in  the  third 
section  is  devoted  to  the  investigation  of  individual  components 
of  the  stress  integral  and  their  dependence  on  the  atmospheric 
response  function  and  the  topography  spectrum  function. 

2. 1  The  Topography  Spectrum  Function 

Very  little  time  was  spent  trying  to  characterize  the 
topography  spectrum  function  in  terms  of  the  features  of  the 
terrain.  We  decided  early  in  the  parameterization  that  the 
spectrum  function  could  be  calculated  and  stored  for  the  grids 
on  which  the  wave  drag  was  to  be  calculated.  The  two-dimensional 
spectrum  of  the  real  topography  is  not  easily  related  to  features 
of  the  terrain. 

The  topography  spectrum  function  is,  by  definition, 

2 

A (k , 1)  =  h*(k,l)  h (k , 1) 

where  X  and  Y  are  the  dimensions  of  the  grid  and  h(k,l)  is 
the  two-dimensional  Fourier  transform  of  the  topography, 
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TABLE  1 

ATMOSPHERES  STUDIED 


ATMOSPHERE  P 

s 

B 

U1 

U3 

V1 

V3 

1 

906.4 

231 

263 

29.3 

7.88 

-0.97 

0.82 

2 

826.0 

234 

262 

33.6 

14.80 

4.57 

-0.69 

3 

906.4 

231 

263 

29.3 

7.96 

-0.97 

0.82 

4 

906.4 

231 

263 

29.3 

8.27 

-0.97 

0.82 

5 

826.0 

234 

262 

33.6 

15.00 

4.57 

-0.69 

6 

826.0 

234 

262 

33.6 

14.60 

4.57 

-0.69 

7 

826.0 

234 

262 

33.6 

15.50 

4.57 

-0.69 

8 

826.0 

234 

262 

33.6 

14.10 

4.57 

-0.69 

9 

826.0 

234 

262 

33.6 

16.30 

4.57 

-0.69 

10 

826.0 

234 

262 

33.6 

13.30 

4.57 

-0.69 

11 

906.4 

231 

263 

33.6 

14.80 

-0.97 

0.82 

12 

906.4 

231 

263 

29.3 

7.88 

4.57 

-0.69 

13 

743.8 

214 

237 

-4.95 

-2.89 

11.47 

5.85 

14 

738.9 

239 

267 

10.95 

5.14 

13.23 

-4.93 

15 

743.8 

239 

267 

— 4 .  r  5 

-2.89 

11.47 

5.85 
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TADLE  2 

TOPOGRAPHIES  INVESTIGATED 


TOPOGRAPHY 

NORTHERN 

LATITUDE 

SOUTHERN 

LATITUDE 

WESTERN 

LONGTIUDE 

EASTERN 

LONGITUDE 

1 

37  °N 

34° 

20  'N 

121°W 

118°  20 ' W 

2 

37  °N 

34° 

2  0  1  N 

118°  20 ' W 

115°  40 'W 

3 

37  °N 

34° 

2  0  1 N 

115°  40 1 W 

113°W 

4 

37  °N 

34  ■* 

2  0  1  N 

113°W 

110°  20 'W 

5 

37  °N 

34° 

20  'N 

110°  20 'W 

107°  40 'W 

6 

37  °N 

34 

20  1 N 

107°  40 1 W 

105°W 

7 

37  °N 

34° 

20  'N 

105  °W 

102°  20 'W 

8 

37  °N 

34° 

2  0 '  N 

102°  20 ' W 

99°  40 'W 

9 

37  °N 

34  ‘ 

20 '  N 

99°  40 '  W 

9  7  °W 

10 

37  °N 

34° 

20  'N 

97  °W 

94°  20 ' W 

h(k,l)  =  ■ 

4 


/ 


Y 

J  h(x,y) 
0 


ei (kx+ly ) 


dxdy  . 


The  function  A(k,l)  uepends  strongly  on  the  amplitude  of  a 
Fourier  component  of  the  topography  and  will  therefore  vary 
over  many  orders  of  magnitude.  The  5'  topography  data  which  we 
have  over  North  America  provides  resolution  of  the  topography 
spectrum  to  a  wavenumber  of  0.3  km  ^  in  the  midlatitudes.  We 
will  consider  the  wavenumber  dependence  of  the  topography  in 
detail  in  our  discussion  of  the  stress  integral  where  the 
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amplitude  of  the  spectrum  functian  must  be  considered  in  the 
light  of  the  dependence  of  the  atmospheric  response  function 
on  wavenumber. 

2.2  The  Atmospheric  Response  Function 

In  the  *>retherton  formulation  of  the  wave  drag  the 
atmospheric  response  function  represents  the  response  of  the 
atmosphere  independent  of  the  forcing  conditions  at  the  lower 
boundary.  If  the  wave  vectors  are  written  in  polar  coordinates 
the  response  of  the  atmosphere  in  the  continuous  regime  has 
been  shown  to  be 


F(k#<J») 


<  dw 
2i  dz 


w*  (z) 


dw 

dz 


w(z) 


w* (0)  w (0) 


(2.1) 


where  w(z)  represents  the  amplitude  of  the  wave  motion.  The 
solution  for  w(z)  is  obtained  by  integrating  Eq.  (1.1)  through 
the  atmosphere  subject  to  the  boundary  conditions  given  in 
Eq.  (1.2).  By  multiplying  the  complex  conjugate  of  Eq.  (1.1) 
by  w  and  subtracting  the  product  of  w*  and  Eq.  (1.1)  it  is 
easily  shown  that  Eq.  (2.1)  is  independent  of  the  z  .  Equa¬ 
tion  (2.1)  can  then  be  simplified  by  evaluating  the  numerator 
at  the  upper  boundary.  The  result, 


F(k,4>) 


2/lT(H)  "  A.,  sgn  [U  (H)  ] 
w‘(0)  +  WI(0) 


(2,2) 


has  been  programmed  into  our  atmospheric  response  function  code. 

From  a  physical  point  of  view  the  atmospheric  response 
function  represents  a  measure  of  the  kinetic  energy  that  must  be 
transferred  from  the  mean  flow  to  the  vertical  at  the  ground 
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level  to  propagate  one  unit  of  kinetic  energy  to  the  level  at 
v'Mich  the  integration  of  the  Scorer  equation  begins.  Reso¬ 
nances  in  the  atmospheric  response  function  will  occur  when  the 
kinetic  energy,  transferred  into  the  vertical  at  the  ground 
level,  required  to  propagate  one  unit  of  energy  into  the 
stratosphere  approaches  zero. 

Equation  (2.2)  illustrates  one  of  the  weaknesses  of  the 
model.  The  atmospheric  response  is  a  function  of  the  height  at 
which  the  integration  of  the  Scorer  equation  begins.  The  de¬ 
pendence  on  the  depth  of  the  atmosphere  appears  explicitly  in 
the  numerator  of  Eq.  (2.2).  of  greater  significance  is  the 
dependence  of  the  denominator  on  the  initial  height,  H,  chosen 
to  begin  the  numerical  integration  of  the  Scorer  equation. 
Typically,  the  stability  is  on  the  order  of  2  x  io-4  sec""2  and 
in  some  layers  of  the  atmosphere,  particularly  near  the  ground, 
the  wind  speed  is  on  the  order  of  1  m/sec.  The  Scorer  param¬ 
eter  is  then  2  x  io  4  m  2  and  large  variations  of  w(0)  will 
then  occur  when  changes  in  z  on  the  order  of  (tt/2)  (1.4  x 
10  m)  or  approximately  100  m  are  made.  Perhaps  because  of 
this  and  for  obvious  reasons  of  computer  economy,  Bretherton  de¬ 
fines  a  critical  layer  at  points  in  the  atmosphere  where  the 
Scorer  parameter,  l2(z)  ,  is  greater  than  2.5  x  io-5  m~2 
His  integration  then  begins  from  this  point  and  no  larger  values 
of  the  Scorer  parameter  occur  in  the  integration. 

To  be  consistent  in  our  studies  we  always  began  the 
integration  at  15  km  above  sea  level.  The  depth  of  the  atmos¬ 
phere  was  then  determined  by  the  average  height  of  the  terrain. 

In  addition  to  those  atmospheres  listed  in  TaLxe  2,  we 
investigated  other  atmospheres  earlier  in  the  contract  period 
to  characterize  the  behavior  of  the  atmospheric  response  func¬ 
tion.  These  atmospheres  were  obtained  by  making  a  linear  inter¬ 
polation  between  the  data  points  taken  from  National  Weather 
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Service  tapes.  Results  of  one  of  these  runs  are  shown  in 
Figures  8  to  13.  The  atmospheric  response  function  was  run  for 
50  valees  of  the  wave  vector  magnitude,  k,  and  for  19  angles 
between  -90°  and  90°  from  the  direction  of  the  zonal  wind.  The 
atmospheric  response  function  is  invariant  under  the  transfor¬ 
mation  from  (k,l)  to  (-k,-l) .  Critical  xayers  for  this  atmos¬ 
phere  occurred  for  angles  around  77°.  Tor  70°  <  4>  <  90°,  the 
Scorer  parameter  is  much  larger  than  the  magnitude  of  the  wave 
vector  and  F(k,l)  is  virtually  independent  of  k  . 

2  2 

When  atmospheric  layers  exist  where  1  (z)  <  k  ,  energy 

will  be  pumped  into  that  mode  which  propagates  to  that  layer 

and  in  phase  with  the  atmospheric  response  at  that  layer.  The 

resonances  which  occur  in  the  continuous  spectrum  in  Figures 

10  to  13  are  associated  with  wave  vectors  such  that  l2(z)  < 

2 

k  for  a  segment  of  the  atmosphere.  These  resonances  will  not 
have  any  effect  on  the  stress  calculated  with  the  topography 
data  that  is  available  to  us.  With  5'  resolution  the  maximum 
resolvable  wavenumber  in  the  topography  spectrum  is  0.3  km-1 
in  the  midlatitudes.  In  the  upper  atmosphere,  where  winds  of 
sufficient  strength  may  be  present  to  make  the  Scorer  parameter 
comparable  to  k2,  the  stability  is  typicallv  4  x  10-^  sec-2. 

The  winds  required  to  give  a  Scorer  parameter  comparable  to  our 
maximum  resolvable  wavenumber  are  on  the  order  of  65  m/sec. 
Without  better  resolution  in  the  atmospheric  variables  it  is 
unlikely  that  the  lapse  rates  or  the  winds  required  to  have  a 
resonance  in  the  atmospheric  response  function  will  be  found. 

In  the  absence  of  these  resonances  the  atmospheric  re¬ 
sponse  is  nearly  independent  of  k  for  a  constant  <J>  .  Over 
the  entire  range  of  ( k ,  <t> )  the  atmospheric  response  changes  by 
less  than  two  orders  of  magnitude  compared  to  the  topography 
spectrum  which  may  vary  over  seven  orders  of  magnitude  or  more. 
These  observations  considerably  simplify  the  treatment  of  the 
stress  integrals  and  make  it  possible  to  neglect  many  of  the 
terms  in  the  integral  of  the  stress. 
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Figure  9.  The  real  and  imaginary  parts  of  the  ground 
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Before  the  beginning  of  the  parameterization  of  the 
wave  drag  it  was  necessary  to  reduce  the  computing  time  re¬ 
quired  for  the  atmospheric  response  function.  ; u  the  start 
of  this  contract  period  the  computing  time  required  to  calcu¬ 
late  the  atmospheric  response  function  for  the  500  wavenumbers , 
required  to  integrate  the  stress  integral,  was  1100  seconds. 

By  a  judicious  optimization  of  the  code  and  a  simplification 
of  the  numerical  methods  involved  in  the  solution  of  the  Scorer 
equation  we  reduced  the  time  required  to  calculate  the  500 
values  of  F(k,<J>)  to  20  seconds  for  our  parameterization  study. 
This  optimization  produced  less  than  a  3%  change  in  the  magni¬ 
tude  of  the  calculated  atmospheric  response  function. 

2.3  The  Calculation  of  the  Stress  Integral 

During  this  contract  period  the  codes  for  calculating 
the  two  components  of  the  sress  were  developed.  The  code  was 
written  to  perform  the  integration  in  horizontal  wave  number 
space;  the  range  and  resolution  of  the  wave  numbers  being  deter¬ 
mined  by  the  resolution  of  the  topography  used  to  calculate  the 
topography  spectrum  function.  The  atmospheric  response  func¬ 
tion  was  calculated  at  each  point  that  the  spectrum  function 
was  calculated.  The  integrals  for  the  two  components  of  the 
stress  are 
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Our  purpose  in  studying  the  stress  integral  was  to  re¬ 
duce  the  number  of  calculations  of  the  atmospheric  response 
function  required  to  calculate  the  stress. 

A  reduction  of  the  time  required  to  solve  the  Scorer 
equation  wil1  not,  by  itself,  simplify  the  wave  drag  calcula¬ 
tion  to  the  extent  necessary  to  meet  the  requirements  of  the 
GCM.  The  actual  computing  time  that  can  be  devoted  to  the 
calculation  of  the  wave  drag  in  the  GCM  will  depend  upon  the 
persistence  of  the  wave  drag  and  the  number  of  cycles  that 
will  not  require  a  recalculation  of  the  drag.  We  have  found 

that  the  wave  drag  is  on  the  order  of  6  to  7  dynes/cm2  on  a 
4°  by  5°  grid  with  the  northwest  corner  located  at  121°W  37°N. 
If  an  upper  limit  of  10  dynes/cm  is  taken  and  we  assume  that 
all  of  the  momentum  is  taken  out  of  the  lower  layer  of  the  GCM 
(typically  3  km  in  depth)  then  the  change  in  the  wind  speed  in 
the  lower  layer  will  be  on  the  order  of  1  m/sec/hour  of  appli¬ 
cation  of  this  stress.  The  meridional  component  of  the  wind 
velocity  is  typically  5  to  10  m/sec  in  the  lower  layer,  hence 
the  wave  drag  should  be  recalculated  at  least  every  hour  in 
the  GCM.  Further  studies  of  the  persistence  of  the  drag  will 
be  performed  on  the  final  version  of  the  wave  drag  subroutine. 

The  running  time  for  older  versions  of  the  GCM  is 
typically  17  min  of  computer  time  per  day  of  calculation.  If 
we  increase  this  running  time  by  5%  for  the  initial  parameteri¬ 
zation  then  2  sec  of  computing  time  will  be  available  for  each 
calculation  if  the  calculation  is  performed  once  in  an  hour. 

Our  calculations  have  shown  that  the  magnitude  of  the  stress 
will  only  be  significant  when  large  mountains  are  present  in  a 
GCM  grid.  We  have  specified  the  grids  that  will  require  a 
calculation  of  the  wave  drag  in  Figure  14.  There  are  on  the 
order  of  400  of  these  grids.  With  2  sec  of  computing  time 
available  for  each  cycle,  we  can  use  5  x  10~2  sec  to  compute 
the  drag  for  each  grid. 
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To  utilize  the  resolution  of  our  5’  topography  data 
to  the  fullest  on  the  GCM  4°  by  3°  grids,  a  64  by  64  grid  is 
required  (the  fast  Fourier  transform  requires  the  number  of 
grids  to  be  powers  of  2).  The  symmetry  properties  of  the  at¬ 
mospheric  response  function  and  the  topography  spectrum  func¬ 
tion,  which  were  discussed  in  the  previous  semiannual  report, 
reduce  the  number  of  grids  required  for  the  calculation  of  the 
stress  integrals  to  32  by  64.  The  simplest  analytic  solution 
of  the  Scorer  equation  is  the  sinusoidal  solution.  The  sine 
function  and  the  cosine  function  require  approximately  30  sec 
of  computing  time  on  the  UCLA  360/91.  In  a  two-level  model 
approximately  10  calculations  of  the  sine  or  cosine  function 
are  required  for  each  wavenumber  component  of  the  stress. 

Hence,  with  the  simplest  analytic  model  fewer  than  16  solu¬ 
tions  of  the  Scorer  equation  can  be  obtained  for  the  calcula¬ 
tion  of  the  two  components  of  the  stress. 

The  first  task  that  was  completed  after  the  develop¬ 
ment  of  the  stress  integral  code  was  to  calculate  the  stress 
integrals  for  a  series  of  topography  grids.  Two  atmospheres 
were  chosen  from  a  run  of  the  Rand  GCM  and  the  stress  was  cal¬ 
culated  over  six  topography  grids.  The  two  atmospheres  are 
listed  as  atmospheres  1  and  2  in  Table  2.  The  first  five 
topography  grids  and  the  last  grid  listed  in  Table  1  were  in¬ 
volved  in  the  calculation.  The  component  of  the  drag  in  the 
direction  of  the  zonal  wind  has  been  illustrated  in  Figure  15 
for  atmosphere  1.  The  results  of  the  drag  calculation  for  the 
other  component  of  the  stress  and  the  two  components  for  the 
other  atmosphere  are  tabulated  in  Table  3.  The  westernmost 
grid  includes  the  Coast  Range,  the  San  Joaquin  Valley  and  parts 
of  the  Sierra  Nevada  range.  It  is  a  grid  with  extreme  varia¬ 
tions  in  height.  The  next  grid  includes  parts  of  the  Sierra 
Nevadas  and  Death  Valley.  The  easternmost  edge  of  T0P05  is 
just  west  of  the  San  Juan  Mountains  and  the  Sangre  de  Cnsto 
Mountains  in  New  Mexico.  The  grids  represented  by  T0P04  and 
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TABLE  3 

THE  RESULTS  OF  THE  STRESS  PARAMETER I Z AT ION 


ATMOSPHERE 

LOCATION 

ZONAL 

STRESS 

MERIDIONAL 

STRESS 

puw/pkA 

pvw/p | i | A 

(puw) (p  £|a) 
(pwP  (pkA) 

Atmosphere  1 

Topography  1 

-12.5 

4.22 

-0.0968 

0.0656 

-1.48 

Topography  2 

-10.6 

4.13 

-0.112 

Topography  3 

-  5.18 

1.38 

-0.0953 

Topography  4 

-  2.95 

1.52 

-0.113 

Topography  5 

-  1.40 

0.596 

-0.11 

Atmosphere  2 

-  0.061 

0.0294 

-0.0795 

0.039 

-2.06 

Topography  1 

13.96 

-5.83 

0.107 

0.0908 

1.18 

Topography  2 

11.62 

-5.59 

0.122 

Topography  3 

5.77 

-3.13 

0.106 

Topography  4 

2.69 

-2.18 

0.103 

Topography  5 

1.39 

-0.92 

0.109 

Atmosphere  3 

Topography  1 

-12.38 

4.15 

-0.096 

0.0646 

-1.48 

Topography  2 

-10.39 

4.17 

-0.109 

Atmosphere  4 

Topography  1 

-11.50 

4.70 

-0.089 

0.0731 

-1.22 

Topography  2 

-  8.54 

3.84 

-0.090 

Atmosphere  5 

Topography  1 

14.90 

-5.28 

0.115 

-0.082 

1.40 

Atmosphere  6 

Topography  1 

13.37 

-5.20 

0.103 

-0.081 

-1.27 

Atmosphere  7 

Topography  1 

15.74 

-5.57 

0.121 

0.087 

1.39 

Atmosphere  8 

Topoqraphy  1 

8.90 

-5.08 

0.0687 

-0.0791 

0.87 

Atmosphere  9 

Topography  1 

15.78 

-6.27 

0.122 

-0.097 

-1.26 

Atmosphere  10 

Topography  1 

-  4.27 

-5.30 

0.033 

-0.082 

-0.402 

jkJ 
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TABLE  3,  continued 


ATMOSPHERE 

LOCATION 

ZONAL 

STRESS 

MERIDIONAL 

STRESS 

puw/pkA 

pvw/p [ 1 1 A 

(puw)  (pUIa) 
(pvw) (pkA) 

Atmosphere  11 
Topography  1 

13.  35 

1.10 

0.103 

0.017 

6.06 

Atmosphere  12 
Topography  1 

-11. 12 

-2.74 

-0.086 

-0.043 

2.0 

Atmosphere  13 
Topography  1 

7.1 

-2.98 

0.055 

-0.046 

1.20 

Atmosphere  14 
Topography  1 

7.73 

-0.48 

0.060 

-7.4*10“ 3 

8.1 

Atmospnere  15 
Topography  1 

2.80 

-1.37 

0.022 

-0.021 

1.05 
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T0P05  are  south  of  the  very  rugged  region  of  the  Rockies  in 

Utah  and  Colorado.  The  last  grid,  TOPOIO,  contains  sections 
of  Arkansas  and  Oklahoma. 

The  early  runs  that  were  made  of  the  copography  spectrum 
function  and  atmospheric  response  function  indicated  that 
F(k,l)  was  nearly  constant  compared  with  the  highly  varying 
spectrum  function.  It  was  found  that  the  topography  spectrum 
function  was  large  for  small  wavenumbers  and  fell  off  rapidly 
as  the  wavenumbers  increased.  If  the  spectrum  function  is 
localized  in  wavenumber  space  then  the  product  U^(0)  F(k,l) 
should  be  nearly  constant  for  the  wavenumbers  which  dominate 
the  stress  integral.  If  this  is  true,  then  dividing  the  stress 
integral  in  Eq.  (2.3)  by 


k  1 
c  c 

2 p  ( 0 )  J  j  kA ( k , 1 )  dkdl 
0  -1 


(2.5) 


and  the  stress  integral  in  Eq.  (2.4)  by 


2p(0>  I  r  HI  A (k ,  1 

0  -1 

c 


)  dkdl 


(2.6) 


should  produce  a  function  which  is  independent  of  the  topography, 
In  Table  3  we  have  tabulated  these  functions  in  the  fourth  and 
fifth  colomns.  The  first  two  atmospheres,  which  were  run  on 
the  five  topographies  indicate  that  the  ratio  of  Eq.  (2.3)  and 
Eq.  (2.5)  has  an  rms  deviation  about  the  mean  of  7%  although  the 
drag  varies  by  an  older  of  magnitude  over  the  five  topographies. 
On  TOPOIO  the  variation  is  more  significant  but  this  can  be 
attributed  to  the  spectrum  function  being  less  localized  be¬ 
cause  there  are  no  significant  features  in  the  topography. 
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We  had  an  error  in  the  integration  routine  for  Eq. 
(2.5)  and  did  not  obtain  much  data  on  this  function  for  the 
meridional  stress.  We  suspect  that  this  function  will  not  be 
as  constant  for  the  meridional  stress  as  for  the  zonal  stress 
because  the  dominant  features  of  the  topography  tend  to  run 
north-south  (with  the  possible  exception  of  the  Himalayas) ; 
therefore,  the  spectrum  function  on  the  1-axis  and  on  off- 
axis  points  is  not  as  large  as  the  values  of  the  spectrum 
function  along  the  k-axis  when  large  mountains  running  north- 
south  are  in  a  grid.  Some  additional  work  is  required  to  de¬ 
cide  how  the  parameterization  of  the  meridional  component  of 
the  stress  will  be  treated. 

In  Table  4  we  have  tabulated  the  wavenumbers  which 
dominate  the  stress  integral.  Most  of  the  zonal  component  of 
the  stress  is  contributed  by  the  wavenumbers  for  which  1=0. 

In  this  case  Un  becomes  the  zonal  wind  component  and  F(k,l) 
is  essentially  constant  if  the  Sccrer  parameter  is  large. 

To  simplify  the  calculation  of  the  stress  it  would  be 

desirable  to  store  the  integral  in  Eq.  (2.4)  for  each  grid  and 

2 

calculate  a  single  value  of  U  F  for  each  atmosphere.  In 

Table  5  we  have  tabulated  the  calculated  value  of  U  F  for 

n 

each  of  the  atmospheres  that  we  have  considered  for  the  wave- 

number  (5.21  x  10  ,  0).  We  have  also  formed  the  product  of 

2 

UnF  with  the  integral  in  Eq.  (2.4)  and  have  compared  it  with 
the  calculated  stress.  For  most  of  the  atmospheres  tabulated 
the  ratio  of  the  stress  to  the  product  ITFpkA  is  between  1.3 
and  2.5.  Hence,  at  least  for  the  meridional  component  of  the 
stress  it  appears  that  the  wave  drag  can  be  calculated  to 
within  approximately  30%  by  calculating  the  atmospheric  re¬ 
sponse  function  for  a  single  wavenumber. 

We  have  also  tabulated,  for  some  of  the  atmospheres, 
the  momentum  flux  which  is  dissipated  at  critical  layers  in 
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the  lower  layer  and  that  which  gets  through  to  the  upper  layer. 
These  results  are  given  in  Table  6.  One  of  the  difficulties 
with  a  two-layer  model  is  deciding  where  the  momentum  which 
propagates  into  the  stratosphere  comes  from.  We  are  making 
the  assumption  that  since  the  momentum  flux  would  not  be 
present  in  the  absence  of  the  mountains  that  the  momentum  must 
be  taken  from  the  lower  layer.  It  is  obvious  that  momentum 
will,  in  fact,  be  taken  from  both  layers  but  it  is  quite  diffi¬ 
cult  to  decide  what  the  contribution  is  from  each  layer. 

3.  PARAMETERIZATION  OF  THE  WAVE  DRAG 

In  this  section  we  summarize  the  progress  that  has  been 
made  toward  delivery  of  a  subroutine  to  the  Rand  Corporation 
for  calculation  of  the  mountain  wave  drag.  We  have  been  de¬ 
layed  by  the  time  required  to  get  acquainted  with  the  DARPA 
network  and  approximately  six  to  eight  man-weeks  of  work  re¬ 
mains  until  the  subroutine  package  is  completed.  Most  of  the 
work  has  been  completed  on  the  parameterization  of  the  zonal 
component  of  the  stress.  During  the  next  phase  of  the  con¬ 
tract  period,  we  will  be  determining  the  number  of  components 
of  the  meridional  stress  that  can  be  calculated  within  the 
limitations  of  the  computing  time  available. 

3 . 1  Work  That  Has  Been  Completed 

In  the  last  section  we  discussed  the  work  that  has  been 
completed  on  the  simplification  of  the  two  components  of  the 
stress.  In  particular,  the  zonal  stress  is  apparently  suscep¬ 
tible  to  great  simplification  over  mountainous  terrain.  We 
propose  to  treat  the  zonal  component  of  the  stress  to  be  inde¬ 
pendent  of  the  meridional  wind  and  to  calculate  a  single  at¬ 
mospheric  response  function  for  the  calculation  of  the  zonal 
stress  integral.  The  results  given  in  Table  5  indicate  that 
thi1"  can  bo  done  (in  particular,  see  atmosphere  13  and  14)  . 
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The  zonal  component  of  the  stress  will  be  determined 
for  each  grid  where  we  anticipate  wave  drag  to  be  important. 
The  atmospheric  response  function  will  be  calculated  analytic- 
aHy  f°r  a  single  wavenumber  on  the  k-axis  in  wavenumber  space 
and  multiplied  by  the  square  of  the  ground  level  zonal  winds. 
This  product  will  be  multiplied  by  the  integral 


kA(k,l)  dkdl  (3.1) 

0  -1 

c 

which  will  be  available  in  core  storage  for  each  grid  to  be 
used  and  a  constant  factor  which  our  results  to  date  indicate 
should  be  approximately  1.9. 

The  flow  diagram  of  the  subroutine  for  the  calculation 
of  the  zonal  stress  is  shown  in  Figure  16.  For  this  initial 
determination  of  the  relative  magnitude  of  the  wave  drag  we  are 
developing  a  two-level  model  of  the  Scorer  parameter.  In  this 
model  constant  values  of  the  Scorer  parameter  are  taken  for 
each  level  of  the  GCM.  The  Scorer  parameter  in  the  upper  level 
of  the  GCM  is  determined  by  the  wind  speed  and  the  temperature 
at  a  =  3/4  .  The  lapse  rate  is  determined  by  the  temperatures 
and  geopotentials  at  the  a  =  1/2  and  o  =  1  levels.  Simi¬ 
larly,  the  Scorer  parameter  can  be  calculated  for  the  lower 
level.  The  assumption  of  a  constant  Scorer  parameter  in  each 
layer  allows  the  simplest  analytic  solution  to  the  Scorer 
equation  to  be  used.  The  initialization  of  the  solution  oc¬ 
curs  at  the  200  mb  level  or  at  a  level  where  a  critical  layer 
occurs.  When  1  (H)  >  k  the  solution  of  the  Scorer  equation 
is  sinusoidal  and  we  take  the  solution  for  outward  propagating 
waves.  In  the  case  where  l2  (H)  <  k2  the  wave  should  be  con¬ 
sidered  a  part  of  the  discrete  spectrum.  However,  in  general. 
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Figure 


16.  Flow  diagram  for  the  calculation  of  the  zonal 
component  of  the  stress. 
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the  trapped  wave  spectrum  will  contribute  little  to  the  total 

drag  and  for  simplicity  of  the  calculation  the  integration  is 

2  2 

started  at  the  level  for  which  1  (H)  >  < 

For  outward  propagating  waves,  the  real  solution  of  the 
Scorer  equation  is 

wR  =  A  cos[(l2(z)  -  <2)z]  (3.2) 


and  the  imaginary  solution  is 


w^  =  B  sin[(l2(z)  -  k2)z] 


(3.3) 


The  boundary  conditions  applied  at  z=H  are  wR=wI=l  ,  cor¬ 
responding  to  one  unit  of  kinetic  energy  arriving  at  this 
level.  The  solution  corresponding  to  this  boundary  condition 
is  found  at  the  a  =  1/2  level  if  z=H  occurs  above  the 
a  =  1/2  level  and  is  then  matched  to  the  corresponding  solu¬ 
tion  in  the  lower  level.  Finally  the  solution  at  the  ground 
level  is  obtained  and  the  ospheric  response  function. 


F (k,l) 


2 (l2 (H)  -  k2)1/2 

2  2  ' 
wR  +  Wl 


(3.4) 


is  evaluated.  Here,  wR  and  w^  are  calculated  at  a  =  1  . 

The  ground  level  winds  are  obtained  from  the  GCM.  The  product 
of  Eqs .  (3.1)  and  (3.4),  the  ground  level  zonal  component  of 
the  wind,  and  the  parameterized  correction  value  give  the  zonal 
stress . 

The  two  level  model  described  above  was  used  by  Palm 
(9) 

and  Foldvik,  with  a  third  level  above  specifying  an  expo¬ 
nentially  decaying  Scorer  parameter.  A  third  analytic  solu- 

(10) 

tion  to  the  Scorer  equation  was  given  by  Vergeiner.  We 
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have  considered  applying  the  Vergeiner  eolation  to  a  two-level 
model  with  a  constant  stability  and  the  wind  profile  that  is 
used  to  extrapolate  the  two-level  winds  in  the  GCM  (essentially 
exponential  in  z) .  The  additional  computing  time  required  to 
determine  the  two  parameters  and  to  match  the  solutions  at  the 
u  -  1/2  level  will  only  bo  justified  if  the  Scorer  parameter 
used  by  Vergeiner  closely  approximates  the  true  atmospheric 
Scorer  parameter.  If  the  actual  Scorer  parameter  profile  is 
denoted  by  lfl(z)  and  the  profile  which  gives  an  analytic 
solution  over  an  atmosphere  of  depth  Az  is  denoted  by  l2(z) 

then  the  condition  for  the  atmospheric  response  to  be  nearly 
the  same  for  the  two  profiles  is 


(2>  ~  k231/2  ds  «  J  ,  (3.5) 


re  Az  z ^  .  When  the  actual  atmosphere  was  formulated 

in  terms  of  the  wind  profile  given  in  the  GCM,  Eq.  (3.5)  could 
not  be  satisfied  with  a  two-level  model  employing  the  Vergeiner 
solution  for  the  two  levels.  Therefore,  the  decision  was  made 
to  use  the  simpler  prescription  for  the  Scorer  parameter. 

In  Table  6  we  have  tabulated  the  vertical  distribution 
of  the  momentum  flux  for  the  atmospheres  that  were  studied.  In 
general,  the  zonal  component  of  the  wind  was  large  in  the  upper 
layer  and  very  little  of  the  momentum  was  dissipated  in  this 
layer.  For  many  of  the  atmospheres,  critical  layers  accounted 
for  the  dissipation  of  very  little  of  the  momentum  flux,  and 
the  majority  of  the  momentum  deficit  in  the  troposphere  propa¬ 
gated  into  the  stratosphere  (most  of  the  flux  listed  under 
puwup ) .  In  the  present  model  the  momentum  flux  which  is  not 
dissipated  in  the  troposphere  is  treated  as  a  deficit  in  the 
lower  layer.  In  reality  there  is  an  exchange  of  momentum  in 


j  [1a(z)  "  <2]1/2  ~  (l2 

s 

21 


50 


SSS-R-74-2331 


both  layers  and  there  should  le  momentum  taken  out  of  both 
layers  of  the  GCM.  We  have  discussed  this  question  in  depth 
with  people  working  in  the  field  without  a  satisfactory  solu¬ 
tion.  This  is  a  fallacy  of  the  linear,  inviscid  model  which 
will  require  an  investigation  of  the  results  of  incorporating 
our  subroutine  into  the  GCM. 

The  output  of  the  wave  drag  subroutine  will  be  two 
forces,  each  to  be  applied  to  one  of  the  two  levels  of  the  GCM. 
The  components  of  the  two  forces  will  act  against  the  meridional 
and  zonal  components  of  the  wind.  We  have  discussed  the  per¬ 
sistence  of  these  forces  in  terms  of  the  effect  these  forces 
hav<_  on  the  wind  components  at  the  two  levels.  A  study  will  be 
made  in  cooperation  with  Mike  Schlesinger  and  Larry  Gates  at 
Rand  to  determine  the  persistence  of  the  drag  in  terms  of  the 
other  forces  acting  on  the  windfield  prior  to  the  incorporation 
of  the  wave  drag  subroutine  into  the  GCM. 

We  have  completed  the  transfer  of  the  topography  tapes 
up  to  UCLA.  These  tapes  have  been  properly  formatted  and  the 
code  has  been  developed  for  calculating  the  topography  spectrum 
function,  and  subsequently,  Eg.  (3.1)  for  each  grid.  We  have 
discovered  gaps  in  the  5'  data  in  those  regions  specified  in 
Figure  14.  Time  will  not  permit  a  thorough  inspection  and  re¬ 
structuring  of  these  tapes;  therefore,  we  will  be  using  only 
the  30'  data  for  all  the  grids  specified  in  Figure  14  on  our 
initial  run  in  the  GCM. 

There  is  an  additional  error  in  the  grids  shown  in 
Figure  14.  We  have  drawn  the  relevant  topography  grids  on  the 
matrix  of  grids  used  in  the  GCM.  These  grids  have  the  velocity 
points  on  the  edges  and  the  scalar  variables  at  the  center  of 
the  grid.  We  will  actually  be  performing  our  calculations  on 
grids  centered  on  velocity  points  and  will  require  four  point 
averages  of  the  scalar  variables  for  each  grid.  Therefore, 
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the  actual  topography  grids  that  will  be  used  are  centered  2° 
to  the  north  and  2° 30'  to  the  west  of  those  grids  shown  in 
Figure  14. 

3 . 2  Work  to  be  Completed  in  this  Contract  Period 

In  order  to  complete  the  parameterization  of  the  wave 
drag  a  series  of  tasks  are  to  be  completed  in  the  next  two 
months.  The  final  task  of  the  contract  will  be  the  delivery  of 
the  subroutine  to  Rand  for  testing  and  incorporation  into  the 
GCM. 

3.2.1  Additional  Runs  of  the  Meridional  Stress  —  The  first 
task  is  being  undertaken  at  this  time.  Additional  runs  of  the 
stress  integral  are  being  made  to  determine  to  what  degree  the 
meridional  stress  can  be  simplified.  V^e  do  not  expect  to  be 
able  to  calculate  a  single  component  of  the  stress  as  we  are 
able  to  do  with  the  zonal  component.  However,  the  simplifica¬ 
tion  that  has  been  found  in  the  zonal  stress  will  provide  addi¬ 
tional  computing  time  for  the  meridional  stress  and  should 
allow  the  calculation  of  the  atmospheric  response  function  for 
ten  to  15  wavenumbers. 

3.2.2  Coding  of  the  Analytic  Solutions  --  The  coding  of  the. 
analytic  solutions  will  follow  a  flow  diagram  similar  to  that 
shown  in  Figure  16.  In  addition,  it  will  contain  a  loop  that 
will  provide  for  a  calculation  of  the  atmospheric  response 
function  for  a  variable  number  of  wavenumbers. 

3.2.3  Determination  of  the  Number  of  Wavenumbers  Used  in  the 
Meridional  Stress  Calculation  —  The  third  task  will  include  a 
series  of  timing  runs  to  determine  how  many  of  the  components 
of  the  Stress  integral  can  be  calculated  in  the  time  allotted. 
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3.2.4  Calculate  the  Spectrum  Function  for  the  Important  Points 
—  The  spectrum  function  will  be  calculated  and  stored  for 
those  wavenumbers  which  will  enter  the  meridional  stress  calcu¬ 
lation.  As  we  have  stated  previously,  the  topography  will  be 
resolved  to  30'  for  this  initial  run  and  only  those  grids 
designated  in  Figure  14  will  be  used.  The  codes  to  calculate 
the  spectrum  function  have  been  developed  and  are  operational. 

3.2  5  Final  Version  of  the  Subroutine  —  After  the  completion 
of  the  above  tasks,  the  final  version  of  the  subroutine  will 
be  written.  For  this  initial  test  on  the  GCM,  the  subroutine 
will  be  designed  to  add  less  than  5%  to  the  running  time  of  the 
global  circulation  model.  The  5%  figure  will  be  based  on  a 
recalculation  of  the  drag  every  hour,  and  the  actual  increase 
in  run  time  will  depend  on  the  persistence  tests  to  be  run  by 
Rand. 

The  code  will  provide  forces  in  mks  units  to  be  applied 
for  the  time  specified  by  Rand. 

3.3  Subsequent  V.Tork  on  the  Wave  Drag  Subroutine 

We  do  not  believe  that  the  wave  drag  subroutine  will 
hcc  been  adequately  studied  at  the  end  of  this  contrect  period. 
Some  additional  work  has  been  discussed  in  the  preceding  sec¬ 
tions  and  we  summarize  here  the  work  that  remains  on  this  pro¬ 
ject.  Much  of  this  work  would  have  been  accomplished  during 
this  contract  period  had  we  not  required  the  changeover  in 
computing  operations. 

3.3.1  Additional  Testing  of  the  Parameterization  —  The 
parameterization  that  is  being  developed  requires  additional 
testing  on  topographies  that  differ  from  those  that  have  been 
used.  In  particular,  we  have  not  adequately  investigated  the 
magnitude  of  the  drag  on  topography  grids  other  than  the  few 
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listed  in  Table  1.  The  lack  of  funds  and  time  has  also  re¬ 
stricted  the  calculations  to  smaller  grids  than  the  4°  by  5° 
grids  used  in  the  GCM.  To  adequately  define  the  parameters 
we  should  perform  additional  calculations  on  the  larger  grids. 

We  would  encourage  the  investigation  of  the  stress 
integrals  over  the  Himalayas  where  the  wave  drag  will  un¬ 
doubtedly  be  large  and  the  components  of  the  topography  spec¬ 
trum  function  may  be  vastly  different  from  those  investigated 
to  date.  In  addition,  the  very  high  winds  that  are  observed 
in  the  GCM  runs  over  the  Himalayas  may  require  a  treatment  of 
the  discrete  wave  spectrum  not  included  in  the  present  param¬ 
eterization. 

3.3.2  Determination  of  Dissipation  Mechanisms  —  Additional 
work  will  be  required  to  determine  at  what  level  momentum  is 
pumped  into  the  waves  and  what  the  means  for  dissipation  of 
the  momentum  carried  by  the  waves  is.  We  have  mentioned  that 
most  of  the  wave  momentum  is  pumped  into  the  very  long  wave¬ 
lengths  which  propagate  into  the  stratosphere.  In  the  present 
parameterization  this  momentum  is  taken  out  of  the  lower  level 
of  the  GCM  and  is  not  put  back  into  the  circulation.  It  may 
be  necessary  to  investigate  the  level  at  which  momentum  is 
removed  by  experimenting  with  the  change  in  the  circulation 
observed  in  the  GCM  as  thi~  momentum  is  removed  from  the  two 
levels  in  varying  proportions. 

Other  dissipation  mechanisms,  not  included  in  the 
present  prescription,  such  as  the  breaking  of  waves  of  large 
amplitude,  should  be  investigated. 

3.3.3  Improvement  of  the  Analytic  Treatment  —  We  have  dis¬ 
cussed  in  some  detail  the  limitations  that  are  placed  on  the 
number  of  calculations  that  can  be  made  to  determine  the  drag. 
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Until  we  have  determined  the  time  required  to  make  the  calcu¬ 
lation  of  the  meridional  component  of  the  wave  drag  it  will  He 
difficult  to  specify  where  improvements  in  the  analytic  treat¬ 
ment  can  be  made. 

4.  RADIATION  IN  THE  EARTH'S  ATMOSPHERE 

The  major  undertakings  for  the  past  contract  period 
were,  first,  to  become  proficient  in  tile  use  of  the  DARPA 
Network  and  the  IBM  360/91  computer  system  at  UCLA;  second, 
to  convert  all  the  atmospheric  radiation  codes  to  run  on  the 
360/91  and  verify  executions  of  the  codes  there  against  cor¬ 
responding  executions  on  the  S3  UNIVAC  1108;  and  third,  to 
make  major  modifications  to  all  the  codes  in  preparation  for 
delivery  to  the  Rand  Climate  Dynamics  Group.  Much  of  this 
work  was  of  a  relatively  routine  nature  but  was  nevertheless 
time-consuming  and  was  made  immeasurably  more  difficult  by  our 
unfamiliarity  with  the  IBM  system,  by  the  lack  of  adequate  con¬ 
sulting  services  at  UCLA,  and  by  the  frequent  system  crashes 
at  UCSD  (our  local  host).  As  a  result,  only  a  few  ATRAD  pro¬ 
duction  runs  have  been  made  on  substantially  new  problems. 

In  the  following,  the  work  which  has  been  accomplished  on 
each  code  to  be  delivered  to  Rand  is  discussed,  followed  by  a 
summary  of  work  which  we  intend  to  complete  by  the  close  of 
the  contract. 

4 . 1  Near-IR  Code 

The  near-IR  code  is  really  a  stand-alone  code  (inde¬ 
pendent  of  ATRAD)  for  calculating  up-flux  and  down-flux  in  an 
aerosol-free,  non-scattering,  non-emitting  atmosphere.  Since 
neither  the  scattering  nor  Planck  terms  enter  the  monochromatic 
radiative  transfer  equation  in  this  case,  a  simple  analytic 
solution  is  possible  for  the  radiative  intensity.  The  down- 
flux  becomes  simply  an  exponential,  and  the  up-flux  an  integral 


55 


SSS-R-74-2331 


over  exponentials  weighted  by  the  directionally-dependent  sur¬ 
face  albedo.  When  these  fluxes  are  integrated  over  a  finite 
spectral  interval,  the  exponentials  are  merely  replaced  by 
transmission  functions.  Hence,  the  calculation  has  only  three 
non-trivial  parts:  the  LOWTRAN '  1  transmission  function 

scheme,  Thekaekara ' s table  of  the  extra-terrestrial  solar 
flux  as  a  function  of  wavelength,  and  an  angular  quadrature  for 

the  up-flux.  Results  from  the  near-IR  code  were  compared  with 

/  7 } 

ATRAD  in  the  last  contract  report  for  the  0.9  —  3.8  y  region 
of  the  spectrum  and  for  two  model  atmospheres  (a  dry  one  ar.d  a 
wet  one) .  It  was  demonstrated  in  that  report  that  there  are 
serious  errors  in  the  Mugge-Moeller  parameterization  of  solar 
absorption  which  is  used  in  the  Mintz-Arakawa  GCM.  According 
to  that  parameterization,  the  fractional  absorption  of  the 
>  0.9  y  solar  flux  due  to  a  vertical  absorber  amount  u*  of 
water  vapor  along  a  slant  path  at  zenith  angle  0  is 


A  (u*sec0) 


0.189 

U-V3., 


(u*sec0 ) 


0.303 


u*  is  calculated  from  a  pressure-scaled  integral  of  the  water 
vapor  density,  and  fQ  is  the  <  0.9  y  fraction  of  the  mean 
annual  extra-terrestrial  solar  flux  SQ  (langleys/  ~,in)  .  It 
was  furthermore  demonstrated  that  no  substantial  improvement 
in  this  parameterization  could  be  obtained  by  using  different 
values  of  the  pressure-scaling  factor  or  of  any  of  the  constants 
in  the  formula.  It  was  suggested  that  an  entirely  new  formula 
be  devised,  one  more  accurately  approximating  the  predictions  of 
both  ATRAD  and  the  near-IR  code. 


In  order  that  the  Rand  Climate  Dynamics  Group  have  the 
facility  to  generate  transmission  data  to  which  to  fit  such  a 
new  formula,  a  copy  of  the  near-IR  code  was  delivered  to  them. 
Some  work  was  also  initiated  here  to  fit  the  solar-flux-weighted 
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transmission  with  piecewise  rational  functions  of  u*  ,  and 
to  that  end  a  rational  function  fitting  code  was  developed. 
Because  of  time  pressures  to  prepare  the  other  radiation  codes 
for  the  360/91,  that  effort  had  to  be  postponed;  however,  we 
are  convinced  that  rational  functions,  being  merely  quotients 
of  polynomials,  and  hence,  computationally  fast  compared  with 
power  laws,  are  the  best  choices  for  fitting  not  only  the  solar 
absorption  function  but  also  many  of  the  other  functions  which 
must  be  parameterized  in  a  GCM.  (Tables  are  also  efficient  in 
terms  of  execution  time,  but  tend  to  swell  core  storage  re¬ 
quirements.)  Our  rational  function  fitting  code  will  be  made 
available  to  Rand  for  use  in  future  parameterization  work. 

4.2  Exponential  Fitting  Tables  Code 

The  exponential  fitting  tables  code  generates  spectral- 
interval-by-spectral-interval  tables  of  exponential-sum  fits 
to  the  transmission,  for  use  by  ATRAD.  The  modifications  which 
have  been  made  to  the  code  are  described  below. 

First,  the  old  scheme  for  combining  terms  with  close 
exponents  in  the  exponential  sum  failed  when  the  code  was  run 
on  the  360/91.  Analysis  of  th ^  problem  revealed  that  the 
double  precision  accuracy  on  the  360/91  was  almost  two  signifi¬ 
cant  digits  less  than  the  UNIVAC  1108  upon  which  the  code  was 
developed,  and  that  the  scheme  (a  secant  method)  was  sensitive 
to  this  decreased  precision.  Rather  than  re-formulate  the 
secant  scheme  to  be  stable,  however,  it  was  replaced  entirely 
with  a  Newton 1 s  method  iteration  which  has  higher-order  con¬ 
vergence  than  the  secant  method  (although  no/  guaranteed  con¬ 
vergence).  The  Newton  scheme  requires  many  fewer  iterations 
than  the  old  secant  method  (generally  3-7  iterations  are  suffi¬ 
cient)  and  in  running  the  exponential  fitting  code  on  a  wide 
variety  of  transmission  data  no  case  has  even  been  encountered 
where  the  new  Newton  scheme  failed  to  converge.  Mathematically, 
the  object  of  the  scheme  is  to  fit  the  sum 
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-k .  u  * 

aie  +  ai+le 


"ki+lu* 


/ 


where  ki  and  ki+1  are  close  (generally  within  5%  of  one  an 
other) ,  with  a  single  term 


The  fit  is  required  to  be  a  least-squares  one,  so  that  the  ob¬ 
ject  is  to  pick  a  and  k  to  minimize  the  residual 

N 

l 

j-1 

The  two  simultaneous  equations  to  be  solved  for  a  and  k  are 
then,  in  the  standard  least-squares  fashion, 

^  P 

31  5  =  0 

=  g  (a  ,k)  =  0 

The  first  of  these  equations  can  be  trivially  solved  for  a 
and  used  to  eliminate  a  in  the  second  equation,  resulting  in 
a  single  (non-linear)  equation  in  k  alone.  This  foregoes, 
however,  the  availability  of  an  excellent  initial  guess  for  a  , 
namely 


a .  e 
i 


-k .  u* 
i  3 


-k. 


+  ai+le 


i+1 


u’ 


-  ae 


-ku* 

3 


a  = 
o 


ai  +  ai+i 


Of  course,  we  also  have  an  excellent  initial  guess  for  k 
which  is 
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ko  ~  I(ki  +  ki+l>  * 

But  the  opportunity  to  use  the  symmetric  initialization  and 
the  advantage  of  having  a  free  to  assume  the  "best"  value  in 
each  iteration  caused  the  two-dimensional  form  of  Newton's 
method  to  be  selected.  The  equations  for  this  iteration  are 

f(an,kn)  +  ff(VV  (an+rV  +  Ie'VV  <kn+l'V  =  0 
9<Vkn>  +  If (an'kn>  (an+l-an>  +  f£(Vkn>  ‘WV  =  0  • 


The  second  group  of  code  changes  concern  the  selection 
of  the  array  of  transmission  function  data  to  be  fitted  in  each 
spectral  interval.  Allowance  has  been  made  for  fitting  sec¬ 
tions  of  a  transmission  curve  (between  absorber-amounts 

Umin  and  umax^  *  is  useful  for  scoping  down  on  parts  of 

trarismission  curves.  The  algorithm  for  selecting  the  array  cf 
transmission  data  T/\v(u*)  then  proceeds  as  follows: 

(1)  If  TAv(umax)  >  Tmax  '  do  not  do  fitting. 


2)  Otherwise,  find  smallest  n  >_  0  such  that 

tAu(u*.  +  0 . 97n (u*  -u*.  ) )  >  T  . 

Av  min  max  mm'  —  mm 

and  let  fitting  interval  be  [u*.  ,U]  ,  where 

mm 

u  =  u*.  +  0.97n(u*  -u*.  )  . 

mm  max  mm' 


(3)  Find  smallest  Ne [Nm .  ,Nm  1  such  that 

min  luaX 

U-u*  . 

T.  (u*  .  )  -  T  (U*  +  min^  <  rn 

Av '  min'  aAv  min  N  -  F  Tdiff 

(if  impossible  to  satisfy  this  criterion,  take  N=N  ) 

max' 
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(4)  The  set  of  transmission  data  is  then 


T4vlu£,in  +  J'1' 


,  N  . 


In  the  discussion  of  this  algorithm,  typical  values  of  the 

parameters  are  given  in  parentheses.  Step  (1)  skips  fitting 

when  the  atmosphere  is  essentially  transparent  (Tmax  =  0.99). 

Step  (2)  assures  that  the  smallest  transmission  value  fitted 

will  be  larger  than  T  .  (T  .  =  0.005).  The  choice  of  T  . 

^  min  min  min 

should  be  guided  by  the  fact  that  the  LOV7TRAN  transmissions 
are  truncated  to  zero  below  0.001.  Step  (3)  attempts  to  keep 
the  difference  between  the  first  and  second  transmission  values 
loss  than  Tdiff  <Tdiff  =  0.01  -  0.02,  Nmin  =  5,  Nmax  =  101). 
Since  the  transmission  curve  becomes  less  and  less  steep  as 
u*  increases,  Tdiff  will  then  (if  the  step  is  successful) 
be  the  largest  difference  between  any  two  transmission  values. 
Ir  the  centers  of  strong  absorption  bands,  step  (3)  is  usually 
unsuccessful  and  in  extreme  cases  the  difference  between  the 
first  two  transmission  values  is  as  large  as  0.5.  This  is  an 
inherent  limitation  in  the  fitting  algorithm,  and  has  two 
causes.  First,  because  of  round-off  error  accumulation,  the 
algorithm  can  only  fit  a  finite  number  of  data  points  (no  more 
than  125-150  with  standard  double  precision  arithmetic) . 

Second,  the  algorithm  requires  data  at  equal  steps  in  u*  in 
order  to  solve  the  linear  least-squares  problem  for  the  coef¬ 
ficients  of  the  exponential  sum  in  a  well-conditioned  way. 

Thus,  there  are  generally  too  many  data  points  for  large  u*  , 
where  the  transmission  varies  slowly,  and  too  few  for  small 
u*  ,  where  it  varies  rapidly.  In  order  to  get  around  the 
equal-steps-of-u*  limitation,  it  would  be  necessary  to  devise 
a  general  least-squares  algorithm  for  the  coefficients  of  an 
arbitrary  sum  of  exponentials.  Unfortunately,  the  latter  is  a 
classical  ill-conditioned  problem  cf  numerical  analysis  (stem¬ 
ming  from  the  fact  that  the  exponentials  form  a  non-orthogonal 
set)  . 


10 


SSS-R-74-2331 


The  third  set  of  code-changes  were  directed  towards  the 
output  formats  of  the  code.  The  code  has  two  forms  of  output: 
first,  the  exponential  fitting  parameters  for  use  by  ATRAD  are 
directed  to  a  (card-image  form)  table;  second,  there  are  prints 
of  the  input  parameters,  the  transmission  function  data  fitted, 
the  corresponding  exponential  fit  and  percent  error,  the  coef¬ 
ficients  and  exponents  of  the  fit,  and  a  summary  of  each  itera¬ 
tion  in  the  fitting  process  if  desired. 

Two  different  forms  for  the  tabular  outt  ut  have  been 
tried.  The  first  dumps  all  the  fitting  parameters  en  masse 
(into  a  namelist  output  format,  although  this  is  obviously  not 
essential) .  The  table  is  then  read  all  at  once  at  the  beginning 
of  ATRAD.  The  advantage  of  this  is  that  only  one  ATRAD  I/O 
swap-out  is  necessary;  the  disadvantage  is  that  tables  of  this 
form  are  not  easily  extended  (since  they  are  read  into  arrays 
or  fixed  dimension)  nor  can  they  be  easily  altered  or  even 
scanned  without  an  interpretive  program.  The  other  babular 
form  dumps  the  fitting  parameters  spectral-interval-by- 
spectral-interval  as  they  are  calculated,  into  a  readable  card- 
image-type  format.  Such  tables  require  several  I/O  requests  by 
ATRAD  for  each  spectral  interval,  but  offer  the  advantage  of 
easy  extensibility,  alterability ,  and  readability.  As  repre¬ 
sentative  CPU  time  requirements  for  ATRAD  in  the  two  cases .  we 
offer  the  following  numbers  from  a  clear-sxy  full-spectrum 
ATRAD  run  for  a  typical  tropical  atmosphere  for  10  sun  angles: 
for  the  "en  masse"  tables,  113  CPU  seconds  (UCLA  IBM  360/91) ; 
for  the  interval  by  interval  tables,  175  CPU  seconds.  There 
is  roughly  a  50%  penalty  in  CPU  time  for  using  the  more 
flexible  table  format.  When  clouds  are  introduced  and  Mie 
tables  must  also  be  read  and  processed,  this  percentage  penalty 
will,  of  course,  be  reduced  substantially. 

The  printed  output  has  been  made  considerably  more 
readable,  particularly  the  prints  of  the  input  variables  con¬ 
trolling  the  calculation.  After  each  input  variable  is  printed 
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out,  its  function  is  now  described  in  detail.  This  feature 
(which  has  also  been  introduced  into  the  other  radiation  codes) 
renders  the  code  much  more  accessible  to  inexperienced  users. 

Two  sets  of  exponential  fitting  tables  were  made  at 
UCLA,  one  for  the  "normal"  ATRAD  spectral  interval  structure 
used  in  all  ATRAD  calculations  to  date,  and  one  for  a  very  re¬ 
fined  spectral  interval  structure  for  calculating  "benchmark 
ATRAD  fluxes  against  which  to  verify  the  fluxes  predicted  with 
the  "normal"  spectral  interval  structure.  The  "normal"  struc¬ 
ture,  using  the  notation  that  nx(An)n2  means  n1  ,nx+An nx+2Ln , 
. . . ,  n2-An,n2,  is  (in  cm  ) : 

60 (60)  600  (20)  800 (40) 1200 (80)  1600(160)  2400  (240)  4  800 

(320) 8000(500) 32000(1000) 35000(1500)48500 

The  refined  structure  is: 

60(20)  2400  (100)  11000  (200)  32000  (500)48500 


The  exponential  fitting  code,  printed  output  from  the 
two  production  runs  ("normal"  and  refined),  and  the  two  fitting 
tables  generated  therewith  are  being  delivered  to  the  Rand 
Climate  Dynamics  Group. 

4.3  ij^+in  Tables  Code 

This  code  makes  tables  of  integrals  over  small  ranges 


Ax  of  size  parameter 

of  the  fundamental 

Mie 

scattering  quanti 

ties,  that  is,  of 

f  a  dx  , 

sea 

/  °ext  dx  ■ 

/ 

(ix+i2)dx 

J 

Ax 

Ax 

Ax 

where 
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M  -  cos  (scattering  angle) 

Vbn'VTn  are  Mie -scattering  functions  which  have 
been  defined  in  a  previous  report 
and  the  notations  are  standard  in  Mie 
scattering  theory. 


It  would,  of  course,  be  desirable  to  construct  defini- 
ables  of  ij+i2  ,  etc.,  but  even  rough  estimates  of  the 
size  of  such  tables  show  that  they  would  bo  far  too  large  to 
be  practical  with  today's  computer  peripheral  storage  and  ac¬ 
cess  times.  The  present  tables  code  makes  tables  of 

the  designated  quantities  only  for  fixed  nr„  and  n.  .  Such 
tables  have  proved  extremely  useful  in  reducing  the'cost  of 
«ie  scattering  computations  for  water  droplet  aerosols  in  the 
visible"  region  of  the  spectrum,  where  n  and  n.  are 
practically  constant.  Even  the  "normal"  ATRAD  spectral  inter¬ 
val  structure  has  44  intervals  between  0.29u  and  0.91u  (in 
order  to  resolve  the  variation  of  Rayleigh  scattering)  and 
the  cost  of  doing  these  44  Mie  calculations  with  re-computatio 

11+12  ■  etc"  «ach  time  would  be  prohibitive.  As  an  ex¬ 
ample  of  the  savings  to  be  obtained  from  using  the  i  +1 

tables,  we  cite  the  following  figures  from  a  test  calculation 
of  the  Mie  tables  code  for  a  typical  Arctic  stratus  cloud  and 
for  three  spectral  intervals  60-120  cm’1,  9500-10000  cm’1 
and  20000-20500  cm'1:  without  the  i  +i  tables  (which 
used  -or  the  20000-20500  cm’1  interval  only),  the  calculation 
required  292  CPU  seconds  on  UCLA(s  IBM  360/91;  with  the  i  +i 
tables,  only  65  CPU  seconds  were  needed.  Estimating  that  sans 
11+i2  tables'  roughly  250  CPU  seconds  were  required  for  the 
20000  20500  cm  interval,  and  taking  this  as  typical  for  the 
0.29-0.91^  region  (which  is  conservative),  we  see  that  the  44 
spectral  intervals  in  this  region  would  require  upwards  of 
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44  ■  ?50  =  11000  CPU  seconds  (3.06  hours)  of  300/91  time.  These 
sort  .  of  estimates  point  up  the  necessity  of  computational  short¬ 
cuts  such  as  the  i-L  +  i2  tables. 

It  should  be  noted  that  the  quantities  tabulated  by  the 

1l  +  i2  code  are  not  asca  '  °ext  '  and  *l  +  i2  themsevles,  but 
integrals  of  these  over  a  small  range  Ax  of  the  size  param- 

«.'^r.  This  alteration  in  the  basic  format  of  the  code  was 
prompted  by  our  inability  to  transfer  the  S3  version  of  the 

tables  to  UCLA  without  a  large  amount  of  conversion,  and 
by  the  realization  that  the  length  of  these  tables  (i  +i 
for  157  scattering  angles  and  for  x=0 . 1  (0 . 1 ) 500  ,  for  a  total 
of  over  785,003  words)  exceeded  the  available  disk  storage  at 
UCLA.  The  tabxes  would  have  had  to  reside  on  tape  and  conse¬ 
quently,  runs  accessing  them  would  in  all  likelihood  have  been 
I/O-bound.  A  reconsideration  of  the  whole  matter  led  to  the 
observation  that  osca  ,  oexfc  ,  and  ix+i2  are  all  used  (in 
the  Mie  tables  code)  .in  integrals  of  the  form 

x 

max 

j  ( i  ^ +i  2  )  n(a)  dx 

x  . 

min 

where  n(a)  is  the  aerosol  size  distribution.  i  +i  *  etc., 
are  very  rapidly  varying  functions  of  x  ,  while  n(a)  =  n(^x) 
is  a  comparatively  slowly  varying  function  of  x  .  Thus,  over 
intervals  Ax  small  with  respect  to  the  scale  of  variation  of 
n  (a)  ,  but  large  with  respect  to  the  scale  of  variation  of 
1^+12  '  we  maY  replace  n(a)  by  some  average  value  n  : 

x+Ax  x+Ax 

j  (li  +  12)  n(a)  dx  =  n(x,Ax)  j  (li  +  12)  dx  • 
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Tiur,  the  full  integra1  may  be  approximated: 


max 


r-  HIM  o 

I  ( i  i +  i  2  )  n 


min 


N 

(a)  dx  =  ^  ^(\in  +  (i-1 ) Ax,  Ax) 
i=l 


/ 


x_. n+iAx 
min 


( i ^  ^  clx  • 


xmin+(i-1)ix 


The  averages  n  are  calculated  in  the  Mie  tables  code  as 


x+Ax 


n  (x.  Ax)  =  ij.  I 


n(a)  dx 


This  ensures  that  n  is  a  representative  value  for  the  inter¬ 
val  Ax  .  Each  value  of  n  is  then  multiplied  by  the  appro¬ 
priate  value  of  f  (i1+i2)  dx  (read  from  the  ix+i2  tables) 

and  the  sum  formed.  The  same  procedure  is  followed  for  a 

and  o  .  sca 

ext 

A  table  of  the  "pre-integrals"  j  (ij+i2)  dx,  etc., 

was  calculated  on  the  360/91  for  n  =  1.335  and  n.  =0  for 

re  im 

use  in  the  0.3-0. 8m  spectral  region  for  water  droplet  aerosols 
A  total  of  157  scattering  angles,  densely  distributed  in  the 
forward  scattering  peak  region,  were  used,  which  were: 

0°(0.1o)2o(0.2o)6o(0.5o)ll°(  l°)20o(2.5o)45o(5°)90o 

plus  the  supplements  of  these  angles.  500  pre-integration 
intervals,  from  (0,1)  to  (499,500)  were  taken,  each  of  width 
x  -  1  .  The  choice  Ax  =  1  was  based  on  the  observation  tha 
the  old  i1+i2  tables  were  made  at  a  0.1  interval  in  x  and 
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that,  therefore,  a  ten-fold 
78,500  words)  was  possible. 
Ax  =  1  will  be  small,  for 
variation  of 


reduction  (from  785,000  words  to 
The  variation  of  n (a)  across 
Ax  =  1  corresponds  to  a  radius 


.  _A  .  X 

A  a=T—  Ax  =  ~- 

2  7T  2  7T 

The  maximum  wavelength  for  which  these  i^t^  tables  might  be 
used  is  A  =  0 . 9y  ,  corresponding  to  Aa  =  0.l4y  ;  and  very 
few  size  distributions  will  vary  much  across  any  0.14y  inter¬ 
val.  Each  pre-integral  was  computed  via  the  trapezoidal  rule 
by  dividing  Ax  =  1  up  into  10  egual  subintervals,  thus  pre¬ 
serving  the  resolution  of  i^t^  t  etc.,  down  to  the  0.1  scale 
demanded  by  Dave.(17)  The  entire  table-making  run  consumed 
1840  CPU  seconds  (0.51  hour)  on  the  IBM  360/91. 

In  order  to  "page  through"  the  i^t^  tables  anc*  scan 
selected  records,  an  interpretive  proaram  mimicing  the  way  the 
Mie  tables  read  the  i1+i2  tables  was  constructed.  This  pro¬ 
gram  verified  the  correctness  of  the  table  described  in  the 
last  paragraph. 

A  number  of  changes  were  made  in  the  code  which  calcu¬ 
lates  o  ,  o  .  ,  and  i,+i0  •  The  most  important  change, 

sea  ext  1  2 

and  one  particularly  suited  to  the  large  core  (4000K  bytes) 
of  the  UCLA  360/91,  is  that  the  functions  ^(y)  and 
are  pre-computed  and  stored  (in  core)  in  very  large  arrays. 
This  very  clearly  leads  to  economies  in  execution  time,  as  op¬ 
posed  to  re-computing  tt^  and  in  for  each  new  value  of  x  . 
However,  the  360/91  essentially  has  free  core  for  jobs  run  on 
the  graveyard  shift;  overall  economies  might  not  be  achieved 
from  pre-computing  and  on  machines  which  charge  for 

core.  (It  might  be  noted  here  that  we  have  also  experimented 
with  keeping  7Tn_Tn  tables  in  peripheral  storage,  and  that 
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this  scheme  is  totally  impractical  given  present-day  access 
times  to  peripheral  storage.) 

Another  code  change  involved  the  elimination  of  complex 
arithmetic.  Complex  arithmetic  would  be  convenient  for  the  Mie 
scattering  calculation,  since  many  of  the  quantities  involved 

(like  a  and  b  )  are  complex.  However,  neither  the  UNIVAC  1108 

n  n 

nor  the  IBM  360/91  FORTRAN  compilers  optimize  complex  opera¬ 
tions  when  one  of  operands  is  real.  For  example,  in  the  divi¬ 
sion  of  a  complex  by  a  real  number,  the  compiler  converts  the 
real  number  to  a  complex  one  (by  adding  a  zero  imaginary  part) 
and  does  a  full  complex  divide,  involving  three  f_oating  adds, 
six  floating  multiplies,  and  two  floating  divides.  If  the 
compiler  recognized  the  denominator  as  real  and  optimized,  only 
two  floating  divides  would  be  necessary.  The  Mie  expressions 
contain  many  instances  where  real  numbers  are  multiplied  or 
divided  with  complex  ones,  and  hence,  these  expressions  are 
unsuitable  for  FORTRAN  complex  arithmetic.  (Considerations 
such  as  these  may  seem  rather  esoteric  to  the  uninitiated,  but 
the  potentially  gigantic  costs  of  Mie  scattering  calculations 
force  one  to  examine  every  statement  in  a  Mie  program  carefully; 
particular  attention  must  be  paid  to  avoiding  re-calculating 
common  expressions,  since  compilers  are  erratic  in  recognizing 
and  storing  such  expressions.) 

Other  less  substantive  code  changes,  particularly  the 
improvement  of  output  formats,  were  also  introduced.  The  en¬ 
tire  code  was  checked  out  against  prior  calculations  and  found 
to  be  in  agreement. 

4 . 4  Mie  Tables  Code 

This  code  calculates  the  integrals  "F  °sca'  °ext' 
and  1^+12  (see  section  4.3)  over  aerosol  size  distribution. 

The  integration  is  done  using  the  trapezoidal  rule  and  is  both 
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printed  and  written  to  unformatted  tables  for  use  by  ATRAD. 

1 1  +  12 '  e^c,»  arG  either  calculated  in-line  or  read  from  tables, 
depending  on  the  wavelength  and  the  values  of  certain  input 
variables . 

A  number  of  changes  needed  to  be  made  in  the  Mie  tables 
code  in  order  to  generalize  it  and  make  it  a  true  production 
code,  for  it  was  the  least  sophisticated  of  the  ATRAC  system. 

It  was  originally  split  off  from  ATRAD  and  retained  many  fea¬ 
tures  which  were  really  irrelevant  to  the  Mie  calculation.  It 
had  a  vertical  zone  structure  mimicing  ATRAD,  and  made  a  Mie 
scattering  calculation  at  each  atmospheric  level  at  which 
aerosol  was  flagged  as  present.  Thus,  a  Mie  table  was  really 
suitable  for  only  one  particular  ATRAD  problem,  and  even  such 
a  simple  change  as  moving  the  aerosol  to  a  different  level  in 
the  ATRAD  calculation  was  non-trivial. 

In  order  to  make  the  Mie  table  as  general  as  possible, 
it  was  necessary  to  allow  it  to  generate  multi-material,  multi¬ 
size-distribution  tables.  All  references  to  level  in  the  at¬ 
mosphere  were  eliminated.  For  each  wavelength,  one  record  is 
now  written  for  each  flagged  material,  and  this  record  contains 
the  polydisperse  Mie  scattering  functions  for  all  of  the  flag¬ 
ged  size  distributions  for  that  material. 

This  general  format  for  the  Mie  tables  allows  a  large 
number  of  different  ATRAD  problems  to  be  run  from  the  same  Mie 
table;  ATRAD  sorts  through  the  data,  taking  only  that  for  the 
materials  and  size  distributions  required  by  its  internal  flans. 
Considerable  computational  economy  is  also  obtained  from  inte- 
gratincr  ij_+i2  '  etc.,  over  many  size  distributions  simultane¬ 
ously,  since  the  most  time-consuming  part  of  the  Mie  calcula¬ 
tion  is  the  evaluation  of  these  mono-disperse  Mie  functions. 

Once  they  are  evaluated,  they  can  be  multiplied  by  any  number 
of  size  distributions  and  added  to  the  trapezoidal  sums  with 
very  little  increase  in  cost  over  just  a  single  size  distribu¬ 
tion  . 
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As  in  the  case  of  the  i , +i2  tables  code  (see  section 
4.3),  a  large-core  version  of  the  Mie  tables  code  was  created 
to  pre-compute  and  save  the  nn's  and  Tn's  needed  for  the  in¬ 
line  calculation  of  i^ij  *  This  version  has  been  checked 
out  against  prior  code  runs  (at  S3)  and  is  working  satisfactorily. 
It  is  approximately  15%  more  efficient  in  terms  of  execution 
rime  than  the  code  version  which  does  not  save  nn  and  . 

A  substantial  amount  of  effort  was  devoted  to  making 
the  output  of  the  Mie  tables  code  self-explanatory  and  complete. 

The  value  and  a  detailed  description  of  each  input  variable  are 
printed  out.  The  available  analytic  size  distributions  (Deir- 
mondjian  modified  gamma,  Gaussian,  log-normal,  and  Junge  power- 
law)  are  edited  with  descriptions  of  how  to  flag  them  and  input 
their  parameters.  The  available  materials  (the  default  set  in¬ 
cludes  water,  ice,  soot,  sea-salt,  and  dust)  are  displayed,  and 
the  tabular  index  of  refraction  of  each  material  actually  used 
is  written  out.  Provision  is  made  for  inputting  the  name  of  a 
material  not  among  the  default  materials  and  its  tabular  index 
of  refraction  as  a  function  of  wavelength,  which  will  be  linearly 
ll\(  jrpolated.  If  i^+i2  tables  are  to  be  used,  the  first  records 
of  those  tables,  containing  identifying  information,  are  read  and 
the  information  edited.  All  of  this  furnishes  a  complete  de¬ 
scription  of  the  problem  setup. 

Code  changes  had  also  to  be  made  to  properly  process 
the  revised  form  of  the  i^t^  tables  (which  now  contain 

I  ( i ^ -*-i 2 )  dx  ,  etc.,  instead  of  simply  i-^+i^  '  etc*^* 

Ax  . 

Furthermore,  the  Henyey-Greenstein  phase  function  option  wnicn 

was  allowed  in  previous  code  versions  has  been  eliminated.  This 

is  because  it  was  found  that  replacing  the  real  phase  function 

by  a  Henyey-Greenstein  approximation,  even  one  having  the  same 

area  under  the  forward  peak,  leads  to  erroneous  flux  predictions 

by  ATRAD.  The  reason  for  this  is  not  so  much  that  the  shape  of 

the  phase  function  is  approximated  (although  the  Henyey-Greenstein 
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phase  function  is  poor  in  this  respect  also) ,  but  that  ATRAD 
truncates  the  forward  peak  of  the  phase  function  if  it  is  high 
enough.  The  truncation  algorithm  often  produces  dramatically 
different  results  for  an  approximating  Henyey-Greenstein  phase 
function  than  for  the  real  phase  function  it  replaces,  because 
the  algorithm  is  tailored  to  the  observed  behavior  of  realistic 
phase  functions.  And  when  a  forward  peak  is  truncated,  the 
truncated  fraction  is  assumed  unscattered  and  the  volume  scatter¬ 
ing  coefficient  accordingly  adjusted  (downwards) .  Thus  the 
approximating  Henyey-Greenstein  phase  function  may  lead  to  a 
scattering  coefficient  quite  different  than  that  calculated  for 
the  real  phase  function,  and  it  is  this  difference  which  primarily 
leads  to  the  large  flux  differences.  (The  Henyey-Greenstein 
option  was  first  introduced  in  order  to  circumvent  the  impossibly 
long  UNIVAC  1108  execution  times  for  the  Mie  tables  code  in  the 
visible  part  of  the  spectrum.  With  the  advent  of  the  i  +i 
tables,  however,  this  problem  was  alleviated  for  water  droplet 
aerosols.  Also,  the  IBM  360/91,  which  is  eight  times  faster 

than  the  UNIVAC  1108 ,  has  brought  all  execution  times  more  within 
reason . ) 

Some  experiments  were  made  at  selected  wavelengths 
(125y,  lp,  and  0.49y)  for  the  Arctic  stratus  size  distribution^ 
to  determine  how  conservative  the  various  parameters  controlling 
the  trapezoidal  integration  over  size  distribution  must  be. 

Two  important  facts  were  discovered: 

(1)  the  observation  of  Dave(17)  that  an  integration 
increment  of  Ax  =  0 . 1  is  sufficient  was  confirmed  (using 
Ax  =  0.05  led  to  only  slight  changes  in  the  results); 

(2)  another  observation  of  Dave's (18),  that  the  inte¬ 
gration  increment  Ax  could  be  increased  to  5Ax  after  a 
large  fraction  of  the  aerosol  particles  (99%)  had  been  inte¬ 
grated  over,  was  proven  wrong;  it  was  necessary  to  use  a 
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larger  traction,  99.9%,  in  order  to  obtain  accurate  results, 
and  with  a  fraction  this  close  to  100%  little  if  any  computa¬ 
tional  savings  are  realized. 


The  second  result  is  not  meant  tc  impugn  the  work  of  Dave,  but 
merely  to  point  up  that  his  method  of  "telescoping"  the  inte¬ 
gration  increment  is  not  general  since  it  did  not  work  for  the 
Arctic  stratus  size  distribution. 


We  offer  a  final  retrospective  observation  on  methods 
of  performing  the  integration  over  size  distribution.  Over 
almost  four  years  we  have  experimented  with  a  number  of  schemes 
to  shorten  the  computation  time  for  this  integration.  These 
schemes  ranged  from  the  sophisticated,  such  as  the  Henyey- 
Greenstein  replacement  (discussed  above)  or  the  Romberg 
extrapolation-to-the-limit  of  multiple  trapezoidal  quadratures, 
down  the  gamut  to  the  simple,  such  as  variants  of  Dave's  tele¬ 
scoping  integration  increment  and  stopping  the  ingration  when 
the  maximum  relative  change  in  the  trapezoidal  sums  in  one 
integration  step  was  small  enough.  All  of  these  schemes  saved 
computer  time,  but  all  of  them  broke  down  and  produced  inac¬ 
curate  results  in  one  situation  or  another.  No  general  scheme 
for  saving  a  substantial  amount  of  computer  time  and  still  pro¬ 
ducing  guaranteed  accurate  results  has  been  found  to  supplant 
the  simple  trapezoidal  integration  from  x  .  to  x  with  a 
step  size 


Ax  =  min 


(0.1, 


X  —  x 

max  min, 
200  ‘ 


4 . 5  ATRAD  —  Code  Improvements 

A  large  number  of  improvements  have  been  made  in  ATRAD 

(7) 

since  the  previous  report  on  this  contract.  Among  the  changes 

which  have  been  introduced  are: 
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U)  a  facility  for  looping  over  multiple  values  of  sun 
angle,  surface  albedo,  surface  temperature,  or 
"top"  temperature; 

(2)  options  for  specular  reflection  and  for  ascribing 
a  temperature  and  emissivity  to  the  top  level  in 
the  mesh  (corresponding,  say,  to  a  cirrus  cloud); 

(3)  greater  efficiency  in  the  Grant-Hunt  algorithm 

/  in  particular  a  cessation  of  calculations 
(when  thermal  emission  is  negligible)  below  the 
level  at  which  sunlight  falls  to  a  ne  ligible 
value  on  account  of  absorption  (applicable  to 
near-IR  and  UV  spectral  intervals) ; 

(4)  the  editing  of  a  detailed  description  of  each 
input  variable,  so  that  ATRAD  is  virtually  self- 
explanatory  in  terms  of  I/O; 

(5)  an  enormous  expansion  of  the  available  options  for 
inputting  atmospheric  vertical  structure; 

(6)  new  logic  for  reading  the  new  generalized  type  of 
Mie  tables,  for  editing  descriptive  information 
from  the  ID  record  of  those  tables,  for  skipping 
through  unwanted  sections  of  the  tables,  and  for 
processing  the  desired  tabular  data  i.ito  the  cor¬ 
rect  atmospheric  levels; 

(7)  allowing  spectral— interval— by-spectral— interval 
flux  edits  at  up  to  four  "slices"  through  the 
vertical  mesh,  thus  alleviating  the  necessity  for 
taking  the  full  flux  edits  at  all  levels  in  order 
to  extract  such  information; 
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(8)  allowing  input  of  tabular  surface  albedos  (diffuse 
and  specular)  as  functions  of  wavelength,  to  be 
linearly  interpolated  to  obtain  albedos  in  each 
spectral  interval;  the  only  alternative  albedo 
option,  a  user-supplied  subroutine  (for  more 
sophisticated  effects  such  as  angular  dependence)  , 
has  been  liberally  supplied  with  comments  to 
assist  the  user;  and, 

(9)  making  edits  more  intelligible  to  the  unfamiliar 
user,  and  supplying  units  information  with  all 
quantities  edited. 

Item  (1),  looping  over  multiple  parameters,  would  be 
trivial  if  it  simply  involved  multiple  executions  of  the  entire 
code.  However,  when  a  parameter  such  as  sun  angle  varies,  many 
of  the  most  time-consuming  parts  of  ATRAD  are  not  affected. 

Thus  what  was  actually  done  was  to  loop  only  over  those  sec¬ 
tions  of  code  involving  the  parameter  in  question.  The  savings 
to  be  achieved  by  doing  this  are  exemplified  by  the  following 
figures:  for  one  sun  angle,  a  full  clear  sky  ATRAD  execution 

(from  60  to  48500  cm  ^)  requires  cbout  40  CPU  seconds  on  the 
IBM  360/91;  for  ten  sun  angles,  this  figure  only  jumps  to  113 
CPU  seconds,  less  than  a  three-fold  increase. 

As  for  the  atmospheric  vertical  structure  (item  (5)), 
the  categories  now  available  are: 

(a)  a  Mintz-Arakawa  2-level  atmosphere,  calculated 
from  temperature,  humidity,  etc.  furnished  for  some  subset  of 
o=0,  1/4,  1/2,  3/4,  1; 

(b)  a  standard  tropical,  midlatitude  (summer  or  winter) , 
or  subarctic  (summer  or  winter)  atmosphere; 
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(c)  a  user-supplied  atmosphere,  for  which  the  sub-options 

are  : 

(i)  to  supply  only  altitude  levels  or  pressure 
levels;  the  missing  variable  will  be  calculated  from 
the  hydrostatic  approximation  taking  account  of 
humidity  and  the  vertical  variation  of  gravity; 

(ii)  to  supply  either  water  vapor  mixing  ratio 
or  water  vapor  density;  the  missing  variable  will  be 
calculated;  and 

(iii)  to  flag  an  analytic  exponentiai-in-alti tude 
or  power-law-in-pressure  mixing  ratio  profile,  with 
parameters  supplied  by  the  user, 

For  both  the  Mintz-Arakawa  and  user-supplied  atmospheres,  there 
are  options  to: 

(a)  flag  an  analytic  ozone  profile,  with  user-input 
parameters ,  and 

(b)  flag  an  exponen tially-decreasing-in-altitude 
aerosol  number  density,  also  with  user-input  parameters. 

For  the  standard  atmospheres ,  the  options  are  for  the  user  to 
input  his  own  altitude  mesh  (which  generates  an  interpolated 
standard  atmosphere),  to  input  his  own  aerosol  structure,  and 
to  scale  the  standard  humidity  and  ozone  profiles  by  constant 
(input)  factors. 

4 . 6  AT  RAD  —  Code  Runs 

ATRAD  is  operational  on  the  UCLA  IBM  360/91  and  has  been 
delivered  to  the  Rand  Climate  Dynamics  Group.  Numerous  ATRAD 
runs  have  been  and  are  being  made  for  Rand  to  demonstrate  ATRAD 
and  its  capabilities. 
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All  ATRAD  runs  to  be  described  below  were  made  for  the 

( 19 ) 

standard  tropical  atmosphere  oi  McClatchey,  et.al.,  for 

the  spectral  region  60-48500  cm  ^  and  for  ten  sun  angles  0sun 

such  that  cos6  =  0.1,  0.2,  ...,  1.0. 
sun 

The  first  pair  of  calculations  which  were  examined  were 
identical  in  all  respects  except  that  one  used  exponential  fit¬ 
ting  tables  generated  some  time  ago  at  S3  on  the  UNIVAC  1108 
and  the  other  used  exponential  fitting  tables  generated  only 
recently  at  UCLA  on  the  IBM  360/91.  The  se  two  tables  were  sub¬ 
stantially  different  because  (a)  the  fitting  is  sensitive  to 
differences  in  machine  precision,  and  (b)  the  method  of  choosing 
transmission  data  to  be  fitted  had  been  changed  vsee  section 

4.2).  Nevertheless,  the  total  fluxes  predicted  in  the  two  runs 

2  2 

never  differed  by  more  than  6  w/m  out  of  1000  w/m  ,  or  0.6%. 

The  heating  rates  between  reasonably  widely-spaced  levels  never 
differed  by  more  than  0.1  degrees/day.  Thus,  the  insensitivity 
of  the  flux  predictions  to  substantial  variations  in  the  expo¬ 
nential  fitting  parameters  is  demonstrated. 

The  second  pair  of  calculations  to  be  compared  were  one 
using  the  "normal"  ATRAD  spectral  interval  structure  and  one 
using  the  much  more  refined  structure  (in  cm  1) 

60  (20) 2400  (100)  11000  (200)  32000(500)48300 

As  discussed  in  section  4.2,  a  second  exponential  fitting  table 
was  made  for  this  refined  structure  on  the  360/91.  Because  of 
the  large  number  of  spectral  intervals  in  the  refined  calcula¬ 
tion  (119  in  the  IR,  4-125u;  84  in  the  near-IR,  0.91-4y;  and 
138  in  the  "visible"  and  UV  0.21-0.91y)  this  calculation  was 
broken  up  into  three  sections,  one  for  the  IR,  one  for  the 
near-IR,  and  one  for  the  "visible"  and  UV.  For  the  IR  part  of 
the  spectrum,  the  fluxes  for  the  "normal"  and  refined  spectral 
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interval  structures  agreed  at  all  vertical  levels  to  within 

±0.3  watts/m'.  For  the  near-IR,  they  agreed  to  within  anywhere 

2  2 
from  ±0.1  watts/m  at  cos0gun  -•  0.1  to  ±1.5  watts/m  at  cos0gun 

=  1.0.  For  the  "visible"  and  UV,  they  agreed  to  within  any- 

2  2 
where  from  ±0.1  watts/m  for  cos0„  =  0.1  to  ±0.3  watts/m  for 

cosBsun  =  1*0*  Thus,  the  "normal"  ATRAD  spectral  interval  struc¬ 
ture  proves  to  be  spectrally  refined  enough  to  satisfy  the  model's 
assumption  that  the  Planck  function,  scattering  and  continuum 
absorption  coefficients,  scattering  phase  function,  and  solar 
flux  do  not  vary  significantly  across  a  spectral  interval  (or 
else,  if  they  do  vary  substantially,  it  is  across  spectral  in¬ 
tervals  which  contribute  negligibly  to  the  flux) . 

2 

The  extreme  variation  of  +1.5  watts/m  in  the  near-IR 
for  a  zenith  sun  deserves  some  further  discussion.  It  is  our 
feeling  that  this  difference,  small  as  it  is  (it  is  only  1/2% 
of  the  flux  in  question) ,  is  too  large  to  be  due  solely  to  the 
spectral  resolution.  We  feel  that  the  bulk  of  the  difference 
is  due  rather  to  the  inadequacy  of  the  exponential  fitting 
algorithm  in  the  strong  near-IR  bands.  As  explained  in 

section  4.2,  that  algorithm  is  simply  unable  to  take  suffici¬ 
ently  many  points  to  resolve  the  transmission  curve  in  strong 
bands.  When  the  spectral  interval  structure  is  more  refined, 
that  problem  is  more  localized  (in  the  band  centers)  and  thus  has 
less  impact  on  the  total  flux.  In  the  future,  a  more  sensible 
ATRAD  spectral  interval  structure  would  take  very  fine  steps 
(perhaps  even  as  small  as  20  cm  i)  across  the  centers  of  these 
strong  bands,  with  much  larger  steps  elsewhere.  In  fact,  this 
measure  had  already  been  taken  in  the  "normal"  ATRAD  spectral 
interval  structure  across  the  15y  CC^  band  (20  cm  steps) 
because  early  in  the  development  of  the  code  it  was  found  that 
unacceptable  flux  errors  were  incurred  when  60  cm  ^  steps  were 
taken . 

The  problem  discussed  in  the  last  paragraph  has  less 
impact  in  the  IR  region,  in  spite  of  the  strong  ^0  bands  there, 
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bocause  in  those  very  bands  where  the  largest  errors  might  be 
incurred,  the  net  fluxes  are  smallest  due  to  the  near-equality 
of  up-  and  down-fluxes. 

A  series  of  ATRAD  runs  have  also  been  made  for  a  wide 
range  of  surface  albedos  (all  the  runs  discussed  above  had  a 
zero  surface  albedo) .  These  will  be  used  to  create  a  table 
showing  the  dependence  of  the  Earth-atmosphere  albedo  on  sun 
angle  and  surface  albedo. 

4 . 7  Future  Work 

Several  tasks  will  be  completed  by  the  close  of  the 
contract.  These  are: 

(1)  the  creation  of  a  Mie  scattering  table  for  water 
for  a  variety  of  aerosol  size  distributions; 

(2)  the  execution  of  ATRAD  for  one  or  two  problems 
using  the  Mie  table  in  (i); 

(3)  delivering  the  following  items  to  Rand  -- 


(a) 

exponential  fitting  code 

(b) 

exponential  fitting  tables  ("normal 
refined) 

"  and 

(c) 

i'l+^2  tables  code 

(d) 

i^+i2  table 

(e) 

Mie  tables  code 

(f) 

Mie  table 

(g) 

rational  function  fitting  code 

(h) 

listings  and,  in  the  case  of  codes, 
executions  of  items  (a)  -  (g); 

sample 

(4)  instructing  Rand  personnel  in  use  of  all  delivered 

codes  ; 
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(5)  complete  work  on  fitting  atmospheric  solar  absorp¬ 
tion  with  rational  functions  of  the  absorber  amount;  and, 

(6)  write  a  summary  of  the  exponential  fitting  al¬ 


gorithm. 
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